 /*	THIS ESTIMATES DEMAND FOR OILS USING HAUSMAN INSTRUMENTS.  IN THIS VARIANT, EACH PRICE IS REGRESSED ON ALL OTHER REGIONS PRICES IN THE FIRST STAGE.  THERE COEFFICIENTS ON EACH REGION'S PRICE IS UNCONSTRAINED*/
set more off
#delimit;
cap log close;
*log using check, replace;
clear all;
set matsize 1000;
*set mem 10000m;
set mem 1200m;
*local path "u:/user11/mw562_rs/my documents/restat/analysis_data";
local path "c:/data/restat_12_9/analysis_data";
local estimator "newey2";
local ivestimator "newey2";
local lags "lag(4)";


local opt "se 3aster bd(2) td(2)";
*sysdir set PERSONAL "u:/user11/mw562_rs/my documents/ado";
di "`path'";

do "c:/data/restat_12_9/programs/csolve_do/csolve";
do "c:/data/restat_12_9/programs/csolve_do/callf";
do "c:/data/restat_12_9/programs/csolve_do/callf_setup";
do "c:/data/restat_12_9/programs/csolve_do/syraids2";


use "`path'/postsyrmaster";
sort REGION BRAND year month day;
bysort REGION BRAND: gen count=_n;


local nbrands1=5;
local nbrands2=4;
local nbrands=9;
local nregions=47;
local list1 "AUNT_JEMIMA HUNGRY_JACK LOG_CABIN MRS_BUTTERWORTH PRIVATE_LABEL";
local list2 "AUNT_JEMIMA_LIGHT JACK_LITE LOG_CABIN_LIGHT MRS_B_LITE";

set seed 10152010;

sort BRAND REGION t;
by BRAND REGION: gen T=_n;

sort T REGION BRAND;
/*separate id for each brand/region combo*/
bysort T: gen id=_n;


/*make list of log price interacted with brand dummies for each segment. exclude interaction with pl because of adding up.  same for instruments*/
foreach y of local list1{;
	forval x=1/4{;
		local reglist1 "`reglist1' lpreg_`x'_`y'";
		local ivreglist1 "`ivreglist1' ivlpreg_`x'_`y'";
	};
};

foreach y of local list2{;
	forval x=1/3{;
		local reglist2 "`reglist2' lplight_`x'_`y'";
		local ivreglist2 "`ivreglist2' ivlplight_`x'_`y'";
	};
};


/*this creates 47 region dummies, one for each of the 5 share equations*/
forval x=2/47{;
	forval y=1/5{;
		gen fe1`x'_`y'=reg`x'*bra1`y';
	};
};
forval x=2/47{;
	forval y=1/4{;
		gen fe2`x'_`y'=reg`x'*bra2`y';
	};
};

forval x=2/47{;
	forval y=1/4{;
		local regionxbrand1 "`regionxbrand1' fe1`x'_`y'";
	};
};

forval x=2/47{;
	forval y=1/3{;
		local regionxbrand2 "`regionxbrand2' fe2`x'_`y'";
	};
};




qui tab month, gen(M);
qui tab year, gen(year);


forval x=1/12{;
	forval y=1/5{;		
		gen M1`x'b`y'=M`x'*bra1`y';
	};
};

forval x=1/12{;
	forval y=1/4{;		
		gen M2`x'b`y'=M`x'*bra2`y';
	};
};

forval x=1/4{;
	gen ivlaidpreg`x'=ivlaidpreg*bra1`x';
};

forval x=1/3{;
	gen ivlaidplight`x'=ivlaidplight*bra2`x';
};



/*this creates month dummies for each of 4 share equations*/
forval x=1/12{;
	forval y=1/4{;
		local Mlist1 "`Mlist1' M1`x'b`y'";
	};
};

forval x=1/12{;
	forval y=1/3{;
		local Mlist2 "`Mlist2' M2`x'b`y'";
	};
};


/*this creates time trends for each of 4 share equations*/
forval x=1/5{;
	gen count1`x'=count*bra1`x';
};

forval x=1/4{;
	gen count2`x'=count*bra2`x';
};



/*UNCONSTRAINED MODEL*/

encode REGION, gen(reg);
tsset id count;
qui compress;

/*HERE ARE THE AIDS OLS ESTIMATES FOR THE REGULAR SEGMENT*/
`estimator' share laidpreg1 laidpreg2 laidpreg3 laidpreg4 `reglist1' `Mlist1' `regionxbrand1' bra11-bra14 count11-count14 if BRAND~="PRIVATE LABEL" & nest==1, nocon `lags';
 

matrix aidB=e(b);
matrix aidVC=e(V);
forval x=1/4{;
	scalar b`x'=_coef[laidpreg`x'];
};

/*recover Private Label share equation parameters through adding up restrictions of AIDS model*/
scalar b5=-b1-b2-b3-b4;
#delimit;
forval x=1/4{;
	foreach y of local list1{;
		scalar g`x'p`y'=_coef[lpreg_`x'_`y'];
	};
};
foreach y of local list1{;
	scalar g5p`y'=-_coef[lpreg_1_`y']-_coef[lpreg_2_`y']-_coef[lpreg_3_`y']-_coef[lpreg_4_`y'];
};
forval x=1/5{;
	scalar g`x'p1=g`x'pAUNT_JEMIMA;
	scalar g`x'p2=g`x'pHUNGRY_JACK;
	scalar g`x'p3=g`x'pLOG_CABIN;
	scalar g`x'p4=g`x'pMRS_BUTTERWORTH;
	scalar g`x'p5=g`x'pPRIVATE_LABEL;
};

forval x=1/`nbrands1'{;
	sum share if bra1`x'==1;
	scalar s1`x'=r(mean);
};




/*HERE ARE THE IV AIDS ESTIMATES FOR THE REGULAR SEGMENT*/
`ivestimator' share `Mlist1' `regionxbrand1' bra11-bra14 count11-count14  (laidpreg1 laidpreg2 laidpreg3 laidpreg4 `reglist1'=ivlaidpreg1 ivlaidpreg2 ivlaidpreg3 ivlaidpreg4 `ivreglist1') if BRAND~="PRIVATE LABEL" & nest==1, nocon `lags';


matrix ivaidB=e(b);
matrix ivaidVC=e(V);
forval x=1/4{;
	scalar ivb`x'=_coef[laidpreg`x'];
};

/*recover Private Label share equation parameters through adding up restrictions of AIDS model*/
scalar ivb5=-ivb1-ivb2-ivb3-ivb4;
#delimit;
forval x=1/4{;
	foreach y of local list1{;
		scalar ivg`x'p`y'=_coef[lpreg_`x'_`y'];
	};
};
foreach y of local list1{;
	scalar ivg5p`y'=-_coef[lpreg_1_`y']-_coef[lpreg_2_`y']-_coef[lpreg_3_`y']-_coef[lpreg_4_`y'];
};

forval x=1/5{;
	scalar ivg`x'p1=ivg`x'pAUNT_JEMIMA;
	scalar ivg`x'p2=ivg`x'pHUNGRY_JACK;
	scalar ivg`x'p3=ivg`x'pLOG_CABIN;
	scalar ivg`x'p4=ivg`x'pMRS_BUTTERWORTH;
	scalar ivg`x'p5=ivg`x'pPRIVATE_LABEL;
};




/*HERE ARE THE OLS AIDS ESTIMATES FOR THE LIGHT SEGMENT*/
`estimator' share laidplight1 laidplight2 laidplight3 `reglist2' `Mlist2' `regionxbrand2' bra21-bra23 count21-count23 if BRAND~="MRS_B_LITE" & nest==2, nocon `lags';
matrix aidB2=e(b);
matrix aidVC2=e(V);
forval x=1/3{;
	scalar b2`x'=_coef[laidplight`x'];
};

/*recover Private Label share equation parameters through adding up restrictions of AIDS model*/
scalar b24=-b21-b22-b23;
#delimit;
forval x=1/3{;
	foreach y of local list2{;
		scalar g2`x'p`y'=_coef[lplight_`x'_`y'];
	};
};
foreach y of local list2{;
	scalar g24p`y'=-_coef[lplight_1_`y']-_coef[lplight_2_`y']-_coef[lplight_3_`y'];
};

forval x=1/4{;
	scalar g2`x'p1=g2`x'pAUNT_JEMIMA_LIGHT;
	scalar g2`x'p2=g2`x'pJACK_LITE;
	scalar g2`x'p3=g2`x'pLOG_CABIN_LIGHT;
	scalar g2`x'p4=g2`x'pMRS_B_LITE;
};

forval x=1/`nbrands2'{;
	sum share if bra2`x'==1;
	scalar s2`x'=r(mean);
};

/*HERE ARE THE AIDS IV ESTIMATES FOR THE LIGHT SEGMENT*/
`ivestimator' share `Mlist2' `regionxbrand2' bra21-bra23 count21-count23 (laidplight1 laidplight2 laidplight3 `reglist2'=ivlaidplight1 ivlaidplight2 ivlaidplight3`ivreglist2' ) /*reg1-reg47*/  if BRAND~="MRS_B_LITE" & nest==2, nocon `lags';

matrix ivaidB2=e(b);
matrix ivaidVC2=e(V);
forval x=1/3{;
	scalar ivb2`x'=_coef[laidplight`x'];
};

/*recover Private Label share equation parameters through adding up restrictions of AIDS model*/
scalar ivb24=-ivb21-ivb22-ivb23;
#delimit;
forval x=1/3{;
	foreach y of local list2{;
		scalar ivg2`x'p`y'=_coef[lplight_`x'_`y'];
	};
};
foreach y of local list2{;
	scalar ivg24p`y'=-_coef[lplight_1_`y']-_coef[lplight_2_`y']-_coef[lplight_3_`y'];
};

forval x=1/4{;
	scalar ivg2`x'p1=ivg2`x'pAUNT_JEMIMA_LIGHT;
	scalar ivg2`x'p2=ivg2`x'pJACK_LITE;
	scalar ivg2`x'p3=ivg2`x'pLOG_CABIN_LIGHT;
	scalar ivg2`x'p4=ivg2`x'pMRS_B_LITE;
};




gen wvol=vol/POP;
forval x=1/5{;
	gen wexpreg`x'=expreg`x'/POP;
};

foreach y of local list1{;
	forval x=1/5{;
		local linplistreg "`linplistreg' preg_`x'_`y'";
	};
};
foreach y of local list1{;
	forval x=1/5{;
		local ivlinplistreg "`ivlinplistreg' ivpreg_`x'_`y'";
	};
};

forval x=2/47{;
	forval y=1/5{;
		local lregionxbrand "`lregionxbrand' fe1`x'_`y'";
	};
};
forval x=1/12{;
	local lMlist1 "`Mlist1' M1`x'b5";
};


/*OLS estimates of linear demand, regular segment*/
`estimator' wvol count11-count15 bra11-bra15 `lMlist1' `lregionxbrand' wexpreg1-wexpreg5 `linplistreg' if nest==1, nocon `lags';  

matrix linB=e(b);
matrix linVC=e(V);
#delimit;
/*p_i_j is coefficient on price j in equation i*/
forval x=1/5{;
	scalar p_`x'_1=_coef[preg_`x'_AUNT_JEMIMA];
	scalar p_`x'_2=_coef[preg_`x'_HUNGRY_JACK];
	scalar p_`x'_3=_coef[preg_`x'_LOG_CABIN];
	scalar p_`x'_4=_coef[preg_`x'_MRS_BUTTERWORTH];
	scalar p_`x'_5=_coef[preg_`x'_PRIVATE_LABEL];
};


forval x=1/5{;
    forval y=1/5{;
        scalar linb`x'p`y'=p_`x'_`y';
    };
};

mat B_reg=J(5,5,0);
forval x=1/5{;
    forval y=1/5{;
        mat B_reg[`x',`y']=linb`x'p`y';
    };
};

mat theta_reg=J(5,1,0);

forval x=1/5{;
	mat theta_reg[`x',1]=_coef[wexpreg`x'];
};	

forval x=1/5{;
	sum price if bra1`x'==1;
	scalar p1`x'=r(mean);
};

forval x=1/`nbrands1'{;
	sum wvol if bra1`x'==1;
	scalar w1`x'=r(mean);
};


`ivestimator' wvol count11-count15 bra11-bra15 `lMlist1' `lregionxbrand' (wexpreg1-wexpreg5 `linplistreg'=ivwexpreg1-ivwexpreg5 `ivlinplistreg') if nest==1, nocon `lags';  


matrix ivlinB=e(b);
matrix ivlinVC=e(V);
#delimit;
forval x=1/5{;
	scalar ivp_`x'_1=_coef[preg_`x'_AUNT_JEMIMA];
	scalar ivp_`x'_2=_coef[preg_`x'_HUNGRY_JACK];
	scalar ivp_`x'_3=_coef[preg_`x'_LOG_CABIN];
	scalar ivp_`x'_4=_coef[preg_`x'_MRS_BUTTERWORTH];
	scalar ivp_`x'_5=_coef[preg_`x'_PRIVATE_LABEL];
};


forval x=1/5{;
    forval y=1/5{;
        scalar ivlinb`x'p`y'=ivp_`x'_`y';
    };
};

mat ivB_reg=J(5,5,0);
forval x=1/5{;
    forval y=1/5{;
        mat ivB_reg[`x',`y']=ivlinb`x'p`y';
    };
};


mat ivtheta_reg=J(5,1,0);

forval x=1/5{;
	mat ivtheta_reg[`x',1]=_coef[wexpreg`x'];
};	

forval x=1/4{;
	gen wexplight`x'=explight`x'/POP;
};


foreach y of local list2{;
	forval x=1/4{;
		local linplistlight "`linplistlight' plight_`x'_`y'";
	};
};

foreach y of local list2{;
	forval x=1/4{;
		local ivlinplistlight "`ivlinplistlight' ivplight_`x'_`y'";
	};
};

forval x=2/47{;
	forval y=1/4{;
		local lregionxbrand2 "`lregionxbrand2' fe2`x'_`y'";
	};
};
forval x=1/12{;
	local lMlist2 "`Mlist2' M2`x'b4";
};



`estimator' wvol count21-count24 bra21-bra24 `lMlist2' `lregionxbrand2' wexplight1-wexplight4 `linplistlight' if nest==2, nocon `lags';  

matrix linB2=e(b);
matrix linVC2=e(V);
#delimit;
forval x=1/4{;
	scalar p2_`x'_1=_coef[plight_`x'_AUNT_JEMIMA_LIGHT];
	scalar p2_`x'_2=_coef[plight_`x'_JACK_LITE];
	scalar p2_`x'_3=_coef[plight_`x'_LOG_CABIN_LIGHT];
	scalar p2_`x'_4=_coef[plight_`x'_MRS_B_LITE];
};


forval x=1/4{;
    forval y=1/4{;
        scalar linb2`x'p`y'=p2_`x'_`y';
    };
};

mat B_light=J(4,4,0);
forval x=1/4{;
    forval y=1/4{;
        mat B_light[`x',`y']=linb2`x'p`y';
    };
};

mat theta_light=J(4,1,0);

forval x=1/4{;
	mat theta_light[`x',1]=_coef[wexplight`x'];
};	


forval x=1/4{;
	sum price if bra2`x'==1;
	scalar p2`x'=r(mean);
};

forval x=1/`nbrands2'{;
	sum wvol if bra2`x'==1;
	scalar w2`x'=r(mean);
};


`ivestimator' wvol bra21-bra24 `lMlist2' `lregionxbrand2'  (wexplight1-wexplight4 `linplistlight'=ivwexplight1-ivwexplight4 `ivlinplistlight') if nest==2, nocon `lags';  


matrix ivlinB2=e(b);
matrix ivlinVC2=e(V);
#delimit;
forval x=1/4{;
	scalar ivp2_`x'_1=_coef[plight_`x'_AUNT_JEMIMA_LIGHT];
	scalar ivp2_`x'_2=_coef[plight_`x'_JACK_LITE];
	scalar ivp2_`x'_3=_coef[plight_`x'_LOG_CABIN_LIGHT];
	scalar ivp2_`x'_4=_coef[plight_`x'_MRS_B_LITE];
};


forval x=1/4{;
    forval y=1/4{;
        scalar ivlinb2`x'p`y'=ivp2_`x'_`y';
    };
};

mat ivB_light=J(4,4,0);
forval x=1/4{;
    forval y=1/4{;
        mat ivB_light[`x',`y']=ivlinb2`x'p`y';
    };
};

mat ivtheta_light=J(4,1,0);

forval x=1/4{;
	mat ivtheta_light[`x',1]=_coef[wexplight`x'];
};	



save temp, replace;


/*Create a dataset where an observation is a region/time/segment market.  Use this to estimate middle level.*/  
drop vol-APRICE merged-lprice lexpreg-ivprice bra11-reg47 ivlaidpreg-ivlaidplight laidpreg1-ivplight_4_MRS_B_LITE ivwexpreg-ivwexplight ivwexpreg1-ivwexplight4 BRAND livprice id ivlaidpreg1-ivlaidplight3 wvol wexpreg1-wexplight4 count11-count24 fe12_1-M212b4;
duplicates drop;


/*Middle segment expenditures*/
gen wexpmid=wexpreg+wexplight;


gen lnvolmid=0;
replace lnvolmid=laidpreg if nest==1;
replace lnvolmid=laidplight if nest==2;
 

tab nest, gen(nest);
gen lp1_1=nest1*indexreg;
gen lp1_2=nest1*indexlight;
gen lp2_1=nest2*indexreg; 
gen lp2_2=nest2*indexlight;

gen p1_1=nest1*lin_indexreg;
gen p1_2=nest1*lin_indexlight;
gen p2_1=nest2*lin_indexreg;
gen p2_2=nest2*lin_indexlight;



gen ivlp1_1=nest1*ivindexreg;
gen ivlp1_2=nest1*ivindexlight;
gen ivlp2_1=nest2*ivindexreg;
gen ivlp2_2=nest2*ivindexlight;

gen ivp1_1=nest1*ivlin_indexreg;
gen ivp1_2=nest1*ivlin_indexlight;
gen ivp2_1=nest2*ivlin_indexreg;
gen ivp2_2=nest2*ivlin_indexlight;

gen lexp1=lexp*nest1;
gen lexp2=lexp*nest2;

gen wexp1=wexp*nest1;
gen wexp2=wexp*nest2;

gen ivlexp1=ivlexp*nest1;
gen ivlexp2=ivlexp*nest2;
gen ivwexp1=ivwexp*nest1;
gen ivwexp2=ivwexp*nest2;


tab month, gen (M);
tab REGION, gen(reg);
forval x=2/47{;
	forval y=1/2{;
		gen fe`y'_`x'=nest`y'*reg`x';
	};
};




forval x=1/12{;
	forval y=1/2{;		
		gen midM`x'b`y'=M`x'*nest`y';
	};
};

#delimit;
forval x=2/47{;
	forval y=1/2{;
		local femid "`femid' fe`y'_`x'";
	};
};

/*this creates month dummies for each of 2 middle level equations*/

forval x=1/12{;
	forval y=1/2{;
		local Mlistmid "`Mlistmid' midM`x'b`y'";
	};
};
gen count1=count*nest1;
gen count2=count*nest2;

sort nest REGION t;
by nest REGION: gen T2=_n;

sort T2 REGION nest;
/*separate id for each brand/region combo*/
bysort T2: gen id=_n;
tsset id T2;

/*For Aids Model*/
`estimator' lnvolmid `femid' `Mlistmid' nest1 nest2 lexp1 lexp2 lp1_1 lp1_2 lp2_1 lp2_2 count1 count2, nocon `lags';
matrix aidMidB=e(b);
matrix aidMidVC=e(V);


sca delta_mid_reg_reg=_coef[lp1_1];
sca delta_mid_reg_light=_coef[lp1_2];
sca delta_mid_light_reg=_coef[lp2_1];
sca delta_mid_light_light=_coef[lp2_2];
sca beta_reg=_coef[lexp1];
sca beta_light=_coef[lexp2];

mat deltamid=[delta_mid_reg_reg, delta_mid_reg_light\delta_mid_light_reg, delta_mid_light_light];
mat betamid=[beta_reg\beta_light];


`ivestimator' lnvolmid `femid' `Mlistmid' nest1 nest2 count1 count2 (lexp1 lexp2 lp1_1 lp1_2 lp2_1 lp2_2=ivlexp1 ivlexp2 ivlp1_1 ivlp1_2 ivlp2_1 ivlp2_2), nocon `lags';


matrix ivaidMidB=e(b);
matrix ivaidMidVC=e(V);

sca ivdelta_mid_reg_reg=_coef[lp1_1];
sca ivdelta_mid_reg_light=_coef[lp1_2];
sca ivdelta_mid_light_reg=_coef[lp2_1];
sca ivdelta_mid_light_light=_coef[lp2_2];
sca ivbeta_reg=_coef[lexp1];
sca ivbeta_light=_coef[lexp2];

mat ivdeltamid=[ivdelta_mid_reg_reg, ivdelta_mid_reg_light\ivdelta_mid_light_reg, ivdelta_mid_light_light];
mat ivbetamid=[ivbeta_reg\ivbeta_light];

/*For linear model*/
`estimator' wexpmid `femid' `Mlistmid' nest1 nest2 wexp1 wexp2 count1 count2 p1_1 p1_2 p2_1 p2_2, nocon `lags';
sca ldelta_mid_reg_reg=_coef[p1_1];
sca ldelta_mid_reg_light=_coef[p1_2];
sca ldelta_mid_light_reg=_coef[p2_1];
sca ldelta_mid_light_light=_coef[p2_2];
sca lbeta_reg=_coef[wexp1];
sca lbeta_light=_coef[wexp2];

mat ldeltamid=[ldelta_mid_reg_reg, ldelta_mid_reg_light\ldelta_mid_light_reg, ldelta_mid_light_light];
mat lbetamid=[lbeta_reg\lbeta_light];


`ivestimator' wexpmid `femid' `Mlistmid' nest1 nest2 count1 count2 (wexp1 wexp2 p1_1 p1_2 p2_1 p2_2=ivwexp1 ivwexp2 ivp1_1 ivp1_2 ivp2_1 ivp2_2), nocon `lags';
sca livdelta_mid_reg_reg=_coef[p1_1];
sca livdelta_mid_reg_light=_coef[p1_2];
sca livdelta_mid_light_reg=_coef[p2_1];
sca livdelta_mid_light_light=_coef[p2_2];
sca livbeta_reg=_coef[wexp1];
sca livbeta_light=_coef[wexp2];

mat livdeltamid=[livdelta_mid_reg_reg, livdelta_mid_reg_light\livdelta_mid_light_reg, livdelta_mid_light_light];
mat livbetamid=[livbeta_reg\livbeta_light];


drop nest lexp-laidplight lin_indexreg lin_indexlight ivlexp ivindexreg ivindexlight ivlin_indexreg ivlin_indexlight ivexp wexpreg-ivwexp2 fe1_2-midM8b2 count1 count2;

duplicates drop;
gen i=0;
forval x=1/47{;
	replace i=`x' if reg`x'==1;
};

sort i t;

bysort i: gen count=_n;
tsset i count;

/*For Aids Model*/
`ivestimator' lqall  M2-M12 reg2-reg47 count (indexall=ivindexall), `lags';
scalar ivdelta_top=_coef[indexall];
matrix iv_b_top=e(b);
matrix iv_vc_top=e(V);
mata: iv_b_top=st_matrix("iv_b_top");
mata: iv_vc_top=st_matrix("iv_vc_top");




`estimator' lqall indexall M2-M12 count reg2-reg47, `lags';
scalar delta_top=_coef[indexall];
matrix ols_b_top=e(b);
matrix ols_vc_top=e(V);
mata: ols_b_top=st_matrix("ols_b_top");
mata: ols_vc_top=st_matrix("ols_vc_top");
mata: mata matsave "c:/data/restat_12_9/analysis_data/post_top_level" ols_b_top ols_vc_top iv_b_top iv_vc_top, replace;


/*For linear model*/
`ivestimator' wexp  M2-M12 reg2-reg47 count (lin_indexall=ivlin_indexall), `lags';
mat ivldelta_top=_coef[lin_indexall];
`estimator' wexp lin_indexall M2-M12 count reg2-reg47, `lags';
mat ldelta_top=_coef[lin_indexall];



#delimit;
use temp, replace;
gen bra1=0;
replace bra1=1 if BRAND=="AUNT_JEMIMA";
gen bra2=0;
replace bra2=1 if BRAND=="HUNGRY_JACK";
gen bra3=0;
replace bra3=1 if BRAND=="LOG_CABIN";
gen bra4=0;
replace bra4=1 if BRAND=="MRS_BUTTERWORTH";
gen bra5=0;
replace bra5=1 if BRAND=="PRIVATE_LABEL";
gen bra6=0;
replace bra6=1 if BRAND=="AUNT_JEMIMA_LIGHT";
gen bra7=0;
replace bra7=1 if BRAND=="JACK_LITE";
gen bra8=0;
replace bra8=1 if BRAND=="LOG_CABIN_LIGHT";
gen bra9=0;
replace bra9=1 if BRAND=="MRS_B_LITE";




gen expmid=expreg+explight;


#delimit;
/*total sales by type*/
/*total syrup sales by type over time*/
sort BRAND REGION nest t;
bysort BRAND REGION nest: egen tsaletype=total(expmid);
/*total syrup sales over time*/
bysort BRAND REGION: egen tsaleall=total(exp);
/*total syrup sales by brand over time*/
bysort BRAND REGION: egen bsale=total(sale);
gen avgsharenest=bsale/tsaletype;
gen avgshare=bsale/tsaleall;


/*w`x'b`y'reg is the average revenue share of producy `y' in region `x' conditional on regular*/
forval x=1/`nregions'{;
    forval y=1/`nbrands1'{;	    
        qui sum avgsharenest if reg`x'==1 & bra1`y'==1;
        scalar w`x'b`y'reg=r(mean);
    };
};

/*w`x'b`y'light is the average revenue share of producy `y' in region `x' conditional on light*/
forval x=1/`nregions'{;
    forval y=1/`nbrands2'{;
        qui sum avgsharenest if reg`x'==1 & bra2`y'==1;
        scalar w`x'b`y'light=r(mean);
    };
};

/*w`x'b`y'all is the average revenue share of producy `y' in region `x'*/
forval x=1/`nregions'{;
    forval y=1/`nbrands'{;
        qui sum avgshare if reg`x'==1 & brand`y'==1;
        scalar w`x'b`y'all=r(mean);
    };
};

/*total sales by segment over all regions and time*/
bysort BRAND nest: egen national_tsale_type=total(expmid);

/*total sales over all regions and time*/
bysort BRAND: egen national_tsale_all=total(exp);


bysort BRAND: egen national_bsale=total(sale);
gen nat_avgshare_nest=national_bsale/national_tsale_type;
gen nat_avgshare_all=national_bsale/national_tsale_all;


forval x=1/`nbrands1'{;
	sum nat_avgshare_nest if bra1`x'==1;
	scalar avgw`x'reg=r(mean);
};

forval x=1/`nbrands2'{;
	sum nat_avgshare_nest if bra2`x'==1;
	scalar avgw`x'light=r(mean);
};

forval x=1/`nbrands1'{;
	sum nat_avgshare_all if bra1`x'==1;
	scalar avgw`x'reg_all=r(mean);
};

forval x=1/`nbrands2'{;
	sum nat_avgshare_all if bra2`x'==1;
	scalar avgw`x'light_all=r(mean);
};



/*HERE WE CONSTRUCT OVER-ALL ELASTICITIES EVALUATED AT AVERAGES TAKEN OVER REGIONS AND TIME*/


/*AIDS estimated by OLS*/
/*regular*/
/*within same nest elasticities*/
forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		scalar ereg`x'`y'=(1+delta_mid_reg_reg)*avgw`y'reg+beta_reg*(1+delta_top)*avgw`y'reg_all+g`x'p`y'/s1`x'+(b`x'/s1`x')*(delta_mid_reg_reg*avgw`y'reg+
		avgw`y'reg_all*(1+delta_top)*beta_reg);
	};
};
forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		scalar elight`x'`y'=(1+delta_mid_light_light)*avgw`y'light+beta_light*(1+delta_top)*avgw`y'light_all+g2`x'p`y'/s2`x'+(b2`x'/s2`x')*(delta_mid_light_light*avgw`y'light+
		avgw`y'light_all*(1+delta_top)*beta_light);
	};
};


forval x=1/`nbrands1'{;
	scalar ereg`x'`x'=-1+ereg`x'`x';
};
forval x=1/`nbrands2'{;
	scalar elight`x'`x'=-1+elight`x'`x';
};
mat regnestelas=J(5,5,0);
mat lightnestelas=J(4,4,0);


forval x=1/5{;
	forval y=1/5{;
	mat regnestelas[`x',`y']=ereg`x'`y';
	};
};

forval x=1/4{;
	forval y=1/4{;
	mat lightnestelas[`x',`y']=elight`x'`y';
	};
};

/*Cross nest elasticities*/
forval x=1/5{;
	forval y=1/4{;
		scalar ereg`x'light`y'=(1+b`x'/s1`x')*(beta_reg*(1+delta_top)*avgw`y'light_all+delta_mid_reg_light*avgw`y'light);
	};
};


mat reg_lightelas=J(5,4,0);

forval x=1/5{;
	forval y=1/4{;
	mat reg_lightelas[`x',`y']=ereg`x'light`y';
	};
};
	

forval x=1/4{;
	forval y=1/5{;
		scalar elight`x'reg`y'=(1+b2`x'/s2`x')*(beta_light*(1+delta_top)*avgw`y'reg_all+delta_mid_light_reg*avgw`y'reg);
	};
};
mat light_regelas=J(4,5,0);

forval x=1/4{;
	forval y=1/5{;
	mat light_regelas[`x',`y']=elight`x'reg`y';
	};
};



mat aidselasticities=J(9,9,0);
mat aidselasticities[1,1]=regnestelas;
mat aidselasticities[6,6]=lightnestelas;
mat aidselasticities[1,6]=reg_lightelas;
mat aidselasticities[6,1]=light_regelas;



forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		scalar ivereg`x'`y'=(1+ivdelta_mid_reg_reg)*avgw`y'reg+ivbeta_reg*(1+ivdelta_top)*avgw`y'reg_all+ivg`x'p`y'/s1`x'+(ivb`x'/s1`x')*(ivdelta_mid_reg_reg*avgw`y'reg+
		avgw`y'reg_all*(1+ivdelta_top)*ivbeta_reg);
	};
};



forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		scalar ivelight`x'`y'=(1+ivdelta_mid_light_light)*avgw`y'light+ivbeta_light*(1+ivdelta_top)*avgw`y'light_all+ivg2`x'p`y'/s2`x'+(ivb2`x'/s2`x')*(ivdelta_mid_light_light*avgw`y'light+
		avgw`y'light_all*(1+ivdelta_top)*ivbeta_light);
	};
};


forval x=1/`nbrands1'{;
	scalar ivereg`x'`x'=-1+ivereg`x'`x';
};
forval x=1/`nbrands2'{;
	scalar ivelight`x'`x'=-1+ivelight`x'`x';
};


mat ivregnestelas=J(5,5,0);
mat ivlightnestelas=J(4,4,0);


forval x=1/5{;
	forval y=1/5{;
	mat ivregnestelas[`x',`y']=ivereg`x'`y';
	};
};

forval x=1/4{;
	forval y=1/4{;
	mat ivlightnestelas[`x',`y']=ivelight`x'`y';
	};
};

/*Cross nest elasticities*/
forval x=1/5{;
	forval y=1/4{;
		scalar ivereg`x'light`y'=(1+ivb`x'/s1`x')*(ivbeta_reg*(1+ivdelta_top)*avgw`y'light_all+ivdelta_mid_reg_light*avgw`y'light);
	};
};
mat ivreg_lightelas=J(5,4,0);

forval x=1/5{;
	forval y=1/4{;
	mat ivreg_lightelas[`x',`y']=ivereg`x'light`y';
	};
};
	

forval x=1/4{;
	forval y=1/5{;
		scalar ivelight`x'reg`y'=(1+ivb2`x'/s2`x')*(ivbeta_light*(1+ivdelta_top)*avgw`y'reg_all+ivdelta_mid_light_reg*avgw`y'reg);
	};
};
mat ivlight_regelas=J(4,5,0);

forval x=1/4{;
	forval y=1/5{;
	mat ivlight_regelas[`x',`y']=ivelight`x'reg`y';
	};
};



mat ivaidselasticities=J(9,9,0);
mat ivaidselasticities[1,1]=ivregnestelas;
mat ivaidselasticities[6,6]=ivlightnestelas;
mat ivaidselasticities[1,6]=ivreg_lightelas;
mat ivaidselasticities[6,1]=ivlight_regelas;



/*Linear Elasticities*/
/*within nest elasticities*/

forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		scalar linereg`x'`y'=(p1`y'/w1`x')*(linb`x'p`y'+theta_reg[`x',1]*(ldeltamid[1,1]*avgw`y'reg+avgw`y'reg_all*lbetamid[1,1]*ldelta_top[1,1]));
	};
};
mat lregelas=J(5,5,0);
forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		mat lregelas[`x',`y']=linereg`x'`y';
	};
};
 
forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		scalar linelight`x'`y'=(p2`y'/w2`x')*(linb2`x'p`y'+theta_light[`x',1]*(ldeltamid[2,2]*avgw`y'light+avgw`y'light_all*lbetamid[2,1]*ldelta_top[1,1]));
	};
};
mat llightelas=J(4,4,0);
forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		mat llightelas[`x',`y']=linelight`x'`y';
	};
};
 
/*across nest elasticities*/
forval x=1/`nbrands1'{;
	forval y=1/`nbrands2'{;
		scalar linereg`x'light`y'=(p2`y'/w1`x')*(theta_reg[`x',1]*(ldeltamid[1,2]*avgw`y'light+lbetamid[1,1]*ldelta_top[1,1]*avgw`y'light_all ));
	};
};
mat lreg_lightelas=J(5,4,0);

forval x=1/`nbrands1'{;
	forval y=1/`nbrands2'{;
		mat lreg_lightelas[`x',`y']=linereg`x'light`y';
	};
};



forval x=1/`nbrands2'{;
	forval y=1/`nbrands1'{;
		scalar linelight`x'reg`y'=(p1`y'/w2`x')*(theta_light[`x',1]*(ldeltamid[2,1]*avgw`y'reg+lbetamid[2,1]*ldelta_top[1,1]*avgw`y'reg_all ));
	};
};

mat llight_regelas=J(4,5,0);

forval x=1/`nbrands2'{;
	forval y=1/`nbrands1'{;
		mat llight_regelas[`x',`y']=linelight`x'reg`y';
	};
};

mat linelasticities=J(9,9,0);
mat linelasticities[1,1]=lregelas;
mat linelasticities[6,6]=llightelas;
mat linelasticities[1,6]=lreg_lightelas;
mat linelasticities[6,1]=llight_regelas;


forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		scalar ivlinereg`x'`y'=(p1`y'/w1`x')*(ivlinb`x'p`y'+ivtheta_reg[`x',1]*(livdeltamid[1,1]*avgw`y'reg+avgw`y'reg_all*livbetamid[1,1]*ivldelta_top[1,1]));
	};
};
mat ivlregelas=J(5,5,0);
forval x=1/`nbrands1'{;
	forval y=1/`nbrands1'{;
		mat ivlregelas[`x',`y']=ivlinereg`x'`y';
	};
};
 
forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		scalar ivlinelight`x'`y'=(p2`y'/w2`x')*(ivlinb2`x'p`y'+ivtheta_light[`x',1]*(livdeltamid[2,2]*avgw`y'light+avgw`y'light_all*livbetamid[2,1]*ivldelta_top[1,1]));
	};
};
mat ivllightelas=J(4,4,0);
forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		mat ivllightelas[`x',`y']=ivlinelight`x'`y';
	};
};
 
/*across nest elasticities*/
forval x=1/`nbrands1'{;
	forval y=1/`nbrands2'{;
		scalar ivlinereg`x'light`y'=(p2`y'/w1`x')*(ivtheta_reg[`x',1]*(livdeltamid[1,2]*avgw`y'light+livbetamid[1,1]*ivldelta_top[1,1]*avgw`y'light_all ));
	};
};
mat ivlreg_lightelas=J(5,4,0);

forval x=1/`nbrands1'{;
	forval y=1/`nbrands2'{;
		mat ivlreg_lightelas[`x',`y']=ivlinereg`x'light`y';
	};
};



forval x=1/`nbrands2'{;
	forval y=1/`nbrands1'{;
		scalar ivlinelight`x'reg`y'=(p1`y'/w2`x')*(ivtheta_light[`x',1]*(livdeltamid[2,1]*avgw`y'reg+livbetamid[2,1]*ivldelta_top[1,1]*avgw`y'reg_all ));
	};
};

mat ivllight_regelas=J(4,5,0);

forval x=1/`nbrands2'{;
	forval y=1/`nbrands1'{;
		mat ivllight_regelas[`x',`y']=ivlinelight`x'reg`y';
	};
};

mat ivlinelasticities=J(9,9,0);
mat ivlinelasticities[1,1]=ivlregelas;
mat ivlinelasticities[6,6]=ivllightelas;
mat ivlinelasticities[1,6]=ivlreg_lightelas;
mat ivlinelasticities[6,1]=ivllight_regelas;


/*s_i: creates average share vector for each region*/
#delimit;
local nregions=47;
local nbrands=9;
forval x=1/`nregions'{;
    mat define sbar`x'=J(`nbrands',1,0);
    forval y=1/`nbrands'{;
        qui sum share if reg`x'==1 & bra`y'==1;
        matrix sbar`x'[`y',1]=_result(3);
	};
};

mat define x_reg=J(47,1,0);
mat define x_light=J(47,1,0);

forval x=1/`nregions'{;
    qui sum expmid if reg`x'==1 & nest==1;
	matrix x_reg[`x',1]=_result(3);
	qui sum expmid if reg`x'==1 & nest==2;
	matrix x_light[`x',1]=_result(3);
};


mat X=[x_reg,x_light];

mat define w_nest=J(47,9,0);

forval x=1/`nregions'{;
	forval y=1/5{;
		mat w_nest[`x',`y']=w`x'b`y'reg;
	};
};
forval x=1/`nregions'{;
	forval y=1/4{;
		mat w_nest[`x',`y'+5]=w`x'b`y'light;
	};
};

mat define w=J(47,9,0);
forval x=1/`nregions'{;
	forval y=1/5{;
		mat w[`x',`y']=w`x'b`y'all;
	};
};
forval x=1/`nregions'{;
	forval y=1/4{;
		mat w[`x',`y'+5]=w`x'b`y'all;
	};
};



/*creates average prices for each region*/
forval x=1/`nregions'{;
    mat define pbar`x'=J(`nbrands',1,0);
    forval y=1/`nbrands'{;
        qui sum price if reg`x'==1 & bra`y'==1;
        matrix pbar`x'[`y',1]=_result(3);
       };
};




/*
mat define tsale=J(`nregions',1,0);
forval x=1/`nregions'{;
    qui sum exp if reg`x'==1;
    matrix tsale[`x',1]=_result(3);
};

gen LTSALE=log(exp);
mat define ltsale=J(`nregions',1,0);
forval x=1/`nregions'{;
    qui sum LTSALE if reg`x'==1;
    matrix ltsale[`x',1]=_result(3);
};


/*creates average log prices for each region*/
forval x=1/`nregions'{;
    mat define lpbar`x'=J(`nbrands',1,0);
    forval y=1/`nbrands'{;
        qui sum lprice if bra`y'==1 & reg`x'==1;
        matrix lpbar`x'[`y',1]=_result(3);
        };
    };
*/	



/*creates average volume for each region*/
forval x=1/`nregions'{;
    mat define vbar`x'=J(`nbrands',1,0);
    forval y=1/`nbrands'{;
        qui sum wvol if bra`y'==1 & reg`x'==1;
        matrix vbar`x'[`y',1]=_result(3);
        };
    };


/*makes pre-merger "ownership" matrices*/
forval i=1/`nregions'{;
    mat define aidomega`i'=J(`nbrands',`nbrands',0);
};

forval i=1/`nregions'{;
    mat define ivaidomega`i'=J(`nbrands',`nbrands',0);
};
forval i=1/`nregions'{;
    mat define X`i'=J(`nbrands',`nbrands',0);
};

/*slope parameters for full solution*/


mat define g_reg=J(`nbrands1',`nbrands1',0);
/*[x,y] element of g_reg is the coefficient on log price of good y in equation x*/
forval x=1/5{;
	forval y=1/5{;
		mat g_reg[`x',`y']=g`x'p`y';
	};
};
mat define g_light=J(`nbrands2',`nbrands2',0);
forval x=1/4{;
	forval y=1/4{;
		mat g_light[`x',`y']=g2`x'p`y';
	};
};

mat define ivg_reg=J(`nbrands1',`nbrands1',0);
forval x=1/5{;
	forval y=1/5{;
		mat ivg_reg[`x',`y']=ivg`x'p`y';
	};
};
mat define ivg_light=J(`nbrands2',`nbrands2',0);
forval x=1/4{;
	forval y=1/4{;
		mat ivg_light[`x',`y']=ivg2`x'p`y';
	};
}; 

mat define b_reg=J(`nbrands1',1,0);
forval x=1/5{;
	mat b_reg[`x',1]=b`x';
};	

mat define b_light=J(`nbrands2',1,0);
forval x=1/4{;
	mat b_light[`x',1]=b2`x';
};	

mat define ivb_reg=J(`nbrands1',1,0);
forval x=1/5{;
	mat ivb_reg[`x',1]=ivb`x';
};	

mat define ivb_light=J(`nbrands2',1,0);
forval x=1/4{;
	mat ivb_light[`x',1]=ivb2`x';
};	
mat delta=[delta_mid_reg_reg, delta_mid_reg_light\ delta_mid_light_reg, delta_mid_light_light];
mat ivdelta=[ivdelta_mid_reg_reg, ivdelta_mid_reg_light\ ivdelta_mid_light_reg, ivdelta_mid_light_light];

mat beta=[beta_reg\beta_light];
mat ivbeta=[ivbeta_reg\ivbeta_light];

mat delta_top=[delta_top];
mat ivdelta_top=[ivdelta_top];
 
forval i=1/`nregions'{;
	forval x=1/`nbrands1'{;
			mat ivaidomega`i'[`x',`x']=(-1+(1+ivdelta_mid_reg_reg)*w_nest[`i',`x']
			+ivbeta_reg*(1+ivdelta_top)*w[`i',`x']
			+ivg`x'p`x'/sbar`i'[`x',1]
			+(ivb`x'/sbar`i'[`x',1])*(ivdelta_mid_reg_reg*w_nest[`i',`x']+w[`i',`x']*(1+ivdelta_top)*ivbeta_reg))*x_reg[`i',1]*sbar`i'[`x',1];
	};
};

forval i=1/`nregions'{;
	forval x=1/`nbrands2'{;
			mat ivaidomega`i'[5+`x',5+`x']=(-1+(1+ivdelta_mid_light_light)*w_nest[`i',5+`x']
			+ivbeta_light*(1+ivdelta_top)*w[`i',5+`x']
			+ivg2`x'p`x'/sbar`i'[5+`x',1]
			+(ivb2`x'/sbar`i'[5+`x',1])*(ivdelta_mid_light_light*w_nest[`i',5+`x']+w[`i',5+`x']*(1+ivdelta_top)*ivbeta_light))*x_light[`i',1]*sbar`i'[`x'+5,1];
	};
};
/*Need to fill out off diagonal of aidomega and ivaidomega matrices*/
forval i=1/`nregions'{;
	mat ivaidomega`i'[6,1]=((1+ivb1/sbar`i'[1,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',6]+ivdelta_mid_reg_light*w_nest[`i',6]))*x_reg[`i',1]*sbar`i'[1,1];
	};
forval i=1/`nregions'{;
	mat ivaidomega`i'[1,6]=((1+ivb21/sbar`i'[6,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',1]+ivdelta_mid_light_reg*w_nest[`i',1]))*x_light[`i',1]*sbar`i'[6,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[7,2]=((1+ivb2/sbar`i'[2,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',7]+ivdelta_mid_reg_light*w_nest[`i',7]))*x_reg[`i',1]*sbar`i'[2,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[2,7]=((1+ivb22/sbar`i'[7,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',2]+ivdelta_mid_light_reg*w_nest[`i',2]))*x_light[`i',1]*sbar`i'[7,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[8,3]=((1+ivb3/sbar`i'[3,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',8]+ivdelta_mid_reg_light*w_nest[`i',8]))*x_reg[`i',1]*sbar`i'[3,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[3,8]=((1+ivb23/sbar`i'[8,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',3]+ivdelta_mid_light_reg*w_nest[`i',3]))*x_light[`i',1]*sbar`i'[8,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[9,4]=((1+ivb4/sbar`i'[4,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',9]+ivdelta_mid_reg_light*w_nest[`i',9]))*x_reg[`i',1]*sbar`i'[4,1];
};
forval i=1/`nregions'{;
	mat ivaidomega`i'[4,9]=((1+ivb24/sbar`i'[9,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',4]+ivdelta_mid_light_reg*w_nest[`i',4]))*x_light[`i',1]*sbar`i'[9,1];
};




forval i=1/`nregions'{;
	forval x=1/`nbrands1'{;
			mat aidomega`i'[`x',`x']=(-1+(1+delta_mid_reg_reg)*w_nest[`i',`x']
			+beta_reg*(1+delta_top)*w[`i',`x']
			+g`x'p`x'/sbar`i'[`x',1]
			+(b`x'/sbar`i'[`x',1])*(delta_mid_reg_reg*w_nest[`i',`x']+w[`i',`x']*(1+delta_top)*beta_reg))*x_reg[`i',1]*sbar`i'[`x',1];
	};
};

forval i=1/`nregions'{;
	forval x=1/`nbrands2'{;
			mat aidomega`i'[5+`x',5+`x']=(-1+(1+delta_mid_light_light)*w_nest[`i',5+`x']
			+beta_light*(1+delta_top)*w[`i',5+`x']
			+g2`x'p`x'/sbar`i'[5+`x',1]
			+(b2`x'/sbar`i'[5+`x',1])*(delta_mid_light_light*w_nest[`i',5+`x']+w[`i',5+`x']*(1+delta_top)*beta_light))*x_light[`i',1]*sbar`i'[`x'+5,1];
	};
};


forval i=1/`nregions'{;
	mat aidomega`i'[6,1]=((1+b1/sbar`i'[1,1])*(beta_reg*(1+delta_top)*w[`i',6]+delta_mid_reg_light*w_nest[`i',6]))*x_reg[`i',1]*sbar`i'[1,1];
	};
forval i=1/`nregions'{;
	mat aidomega`i'[1,6]=((1+b21/sbar`i'[6,1])*(beta_light*(1+delta_top)*w[`i',1]+delta_mid_light_reg*w_nest[`i',1]))*x_light[`i',1]*sbar`i'[6,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[7,2]=((1+b2/sbar`i'[2,1])*(beta_reg*(1+delta_top)*w[`i',7]+delta_mid_reg_light*w_nest[`i',7]))*x_reg[`i',1]*sbar`i'[2,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[2,7]=((1+b22/sbar`i'[7,1])*(beta_light*(1+delta_top)*w[`i',2]+delta_mid_light_reg*w_nest[`i',2]))*x_light[`i',1]*sbar`i'[7,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[8,3]=((1+b3/sbar`i'[3,1])*(beta_reg*(1+delta_top)*w[`i',8]+delta_mid_reg_light*w_nest[`i',8]))*x_reg[`i',1]*sbar`i'[3,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[3,8]=((1+b23/sbar`i'[8,1])*(beta_light*(1+delta_top)*w[`i',3]+delta_mid_light_reg*w_nest[`i',3]))*x_light[`i',1]*sbar`i'[8,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[9,4]=((1+b4/sbar`i'[4,1])*(beta_reg*(1+delta_top)*w[`i',9]+delta_mid_reg_light*w_nest[`i',9]))*x_reg[`i',1]*sbar`i'[4,1];
};
forval i=1/`nregions'{;
	mat aidomega`i'[4,9]=((1+b24/sbar`i'[9,1])*(beta_light*(1+delta_top)*w[`i',4]+delta_mid_light_reg*w_nest[`i',4]))*x_light[`i',1]*sbar`i'[9,1];
};

forval i=1/`nregions'{;
    forval x=1/`nbrands1'{;
        mat X`i'[`x',`x']=x_reg[`i',1];
    };
};

forval i=1/`nregions'{;
    forval x=1/`nbrands2'{;
        mat X`i'[`x'+5,`x'+5]=x_light[`i',1];
    };
};



/*
forval x=1/`nbrands2'{;
	forval y=1/`nbrands2'{;
		scalar ivelight`x'`y'=((1+ivdelta_mid_light_light)*avgw`y'light+ivbeta_light*(1+ivdelta_top)*avgw`y'light_all+ivg2`x'p`y'/s2`x'+(ivb2`x'/s2`x')*(ivdelta_mid_light_light*avgw`y'light+
		avgw`y'light_all*(1+ivdelta_top)*ivbeta_light));
	};
};


forval x=1/`nbrands1'{;
	scalar ivereg`x'`x'=-1+ivereg`x'`x';
};
forval x=1/`nbrands2'{;
	scalar ivelight`x'`x'=-1+ivelight`x'`x';
};


forval x=1/`nbrands'{;
    forval i=1/`nregions'{;
        mat aidomega`i'[`x',`x']=-1+((g`x'p`x'-w`i'b`x'*b`x')/sbar`i'[`x',1])+(1+`e')*sbar`i'[`x',1];
    };
};



forval i=1/`nbrands'{;
    forval j=1/`nbrands'{;
        mat g[`i',`j']=g`i'p`j';
    };
};

forval i=1/`nbrands'{;
    forval j=1/`nbrands'{;
        mat ivg[`i',`j']=ivg`i'p`j';
    };
};


forval i=1/`nbrands'{;
	mat b[`i',1]=b`i';
};

forval i=1/`nbrands'{;
	mat ivb[`i',1]=ivb`i';
};

forval i=1/`nregions'{;
    forval j=1/`nbrands'{;
        mat w[`i',`j']=w`i'b`j';
    };
};
    
forval x=1/`nbrands'{;
    forval i=1/`nregions'{;
        mat aidomega`i'[`x',`x']=-1+((g`x'p`x'-w`i'b`x'*b`x')/sbar`i'[`x',1])+(1+`e')*sbar`i'[`x',1];
    };
};

forval x=1/`nbrands'{;
    forval i=1/`nregions'{;
        mat ivaidomega`i'[`x',`x']=-1+((ivg`x'p`x'-w`i'b`x'*ivb`x')/sbar`i'[`x',1])+(1+`e')*sbar`i'[`x',1];
    };
};



/*computes pre-merger markups for each of the ten regions*/
forval i=1/`nregions'{;
    mat ds`i'=diag(sbar`i');
};
*/



/*aidmu has p-mc, not lerner index*/
forval i=1/`nregions'{;
    mat aidmu`i'=J(`nbrands',1,0);
    mat aidmu`i'=-inv(aidomega`i')*X`i'*sbar`i';
    mat ivaidmu`i'=J(`nbrands',1,0);
    mat ivaidmu`i'=-inv(ivaidomega`i')*X`i'*sbar`i';
};



/*mark-up matrix has 10 columns, which correspond to regional markets, and 7 rows corresponding to brands.  
aidmu are p-mc's with aid system, lmu are p=mc's with linear system, llmu are lerner indices with log-linear system
this just puts them in matrix form as described above*/
mat aidmu=[aidmu1];
mat ivaidmu=[ivaidmu1];
mat pbar=[pbar1];
/*mat lpbar=[lpbar1];*/
mat sbar=[sbar1];
forval x=2/`nregions'{;
    mat aidmu=[aidmu,aidmu`x'];
    mat ivaidmu=[ivaidmu,ivaidmu`x'];
    mat pbar=[pbar,pbar`x'];
    mat sbar=[sbar,sbar`x'];
};

mat pbar=pbar';
/*mat lpbar=lpbar';*/
mat sbar=sbar';
mat aidmu=aidmu';
mat ivaidmu=ivaidmu';

/*get costs by solving mark-ups.  because element-by-element operations are a pain in stata, i do this in mata*/
/*input markup matrices, average prices and shares, and post-merger ownership matrices into stata*/
mata;
	b_reg=st_matrix("b_reg");
	b_light=st_matrix("b_light");
	ivb_reg=st_matrix("ivb_reg");
	ivb_light=st_matrix("ivb_light");
	
	g_reg=st_matrix("g_reg");
	g_light=st_matrix("g_light");
	ivg_reg=st_matrix("ivg_reg");
	ivg_light=st_matrix("ivg_light");
	w=st_matrix("w");
	w_nest=st_matrix("w_nest");
	X=st_matrix("X");
	delta=st_matrix("delta");
	ivdelta=st_matrix("ivdelta");
	
	beta=st_matrix("beta");
	ivbeta=st_matrix("ivbeta");
	
	delta_top=st_matrix("delta_top");
	ivdelta_top=st_matrix("ivdelta_top");
	
    aidmu=st_matrix("aidmu");
    /*ivb=st_matrix("ivb");
    ivg=st_matrix("ivg");*/
    ivaidmu=st_matrix("ivaidmu");

    pbar=st_matrix("pbar");
    sbar=st_matrix("sbar");
    aidcost=pbar:*(J(`nregions',`nbrands',1)-aidmu);
    ivaidcost=pbar:*(J(`nregions',`nbrands',1)-ivaidmu);
    st_matrix("aidcost", aidcost);
    st_matrix("ivaidcost", ivaidcost);
end;



/*defines post-merger ownership matrix of elasticities*/
forval i=1/`nregions'{;
    mat define taidomega`i'=J(`nbrands',`nbrands',0);
};
forval i=1/`nregions'{;
    mat define ivtaidomega`i'=J(`nbrands',`nbrands',0);
};




forval i=1/`nregions'{;
	mat ivtaidomega`i'[8,4]=((1+ivb4/sbar`i'[4,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',8]+ivdelta_mid_reg_light*w_nest[`i',8]))*x_reg[`i',1]*sbar`i'[4,1];
	};
forval i=1/`nregions'{;
	mat ivtaidomega`i'[4,8]=((1+ivb23/sbar`i'[8,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',4]+ivdelta_mid_light_reg*w_nest[`i',4]))*x_light[`i',1]*sbar`i'[8,1];
};

forval i=1/`nregions'{;
	mat ivtaidomega`i'[9,3]=((1+ivb3/sbar`i'[3,1])*(ivbeta_reg*(1+ivdelta_top)*w[`i',9]+ivdelta_mid_reg_light*w_nest[`i',9]))*x_reg[`i',1]*sbar`i'[3,1];
};
forval i=1/`nregions'{;
	mat ivtaidomega`i'[3,9]=((1+ivb24/sbar`i'[9,1])*(ivbeta_light*(1+ivdelta_top)*w[`i',9]+ivdelta_mid_light_reg*w_nest[`i',3]))*x_light[`i',1]*sbar`i'[9,1];
};



forval i=1/`nregions'{;            
			mat ivtaidomega`i'[4,3]=((1+ivdelta_mid_reg_reg)*w_nest[`i',4]
			+ivbeta_reg*(1+ivdelta_top)*w[`i',4]
			+ivg3p4/sbar`i'[3,1]
			+(ivb3/sbar`i'[3,1])*(ivdelta_mid_reg_reg*w_nest[`i',4]+w[`i',4]*(1+ivdelta_top)*ivbeta_reg))*x_reg[`i',1]*sbar`i'[3,1];
	};
	

forval i=1/`nregions'{;
			mat ivtaidomega`i'[3,4]=((1+ivdelta_mid_reg_reg)*w_nest[`i',3]
			+ivbeta_reg*(1+ivdelta_top)*w[`i',3]
			+ivg4p3/sbar`i'[4,1]
			+(ivb4/sbar`i'[4,1])*(ivdelta_mid_reg_reg*w_nest[`i',3]+w[`i',3]*(1+ivdelta_top)*ivbeta_reg))*x_reg[`i',1]*sbar`i'[4,1];
};


forval i=1/`nregions'{;
			mat ivtaidomega`i'[9,8]=((1+ivdelta_mid_light_light)*w_nest[`i',9]
			+ivbeta_light*(1+ivdelta_top)*w[`i',9]
			+ivg23p4/sbar`i'[8,1]
			+(ivb23/sbar`i'[8,1])*(ivdelta_mid_light_light*w_nest[`i',9]+w[`i',9]*(1+ivdelta_top)*ivbeta_light))*x_light[`i',1]*sbar`i'[8,1];
};

forval i=1/`nregions'{;
	mat ivtaidomega`i'[8,9]=((1+ivdelta_mid_light_light)*w_nest[`i',8]
			+ivbeta_light*(1+ivdelta_top)*w[`i',8]
			+ivg24p3/sbar`i'[9,1]
			+(ivb24/sbar`i'[9,1])*(ivdelta_mid_light_light*w_nest[`i',8]+w[`i',8]*(1+ivdelta_top)*ivbeta_light))*x_light[`i',1]*sbar`i'[9,1];
};


forval i=1/`nregions'{;
	mat taidomega`i'[8,4]=((1+b4/sbar`i'[4,1])*(beta_reg*(1+delta_top)*w[`i',8]+delta_mid_reg_light*w_nest[`i',8]))*x_reg[`i',1]*sbar`i'[4,1];
	};


forval i=1/`nregions'{;
	mat taidomega`i'[4,8]=((1+b23/sbar`i'[8,1])*(beta_light*(1+delta_top)*w[`i',4]+delta_mid_light_reg*w_nest[`i',4]))*x_light[`i',1]*sbar`i'[8,1];
};

forval i=1/`nregions'{;
	mat taidomega`i'[9,3]=((1+b3/sbar`i'[3,1])*(beta_reg*(1+delta_top)*w[`i',9]+delta_mid_reg_light*w_nest[`i',9]))*x_reg[`i',1]*sbar`i'[3,1];
};


forval i=1/`nregions'{;
	mat taidomega`i'[3,9]=((1+b24/sbar`i'[9,1])*(beta_light*(1+delta_top)*w[`i',3]+delta_mid_light_reg*w_nest[`i',3]))*x_light[`i',1]*sbar`i'[9,1];
};



forval i=1/`nregions'{;            
			mat taidomega`i'[4,3]=((1+delta_mid_reg_reg)*w_nest[`i',4]
			+beta_reg*(1+delta_top)*w[`i',4]
			+g3p4/sbar`i'[3,1]
			+(b3/sbar`i'[3,1])*(delta_mid_reg_reg*w_nest[`i',4]+w[`i',4]*(1+delta_top)*beta_reg))*x_reg[`i',1]*sbar`i'[3,1];
	};

	
forval i=1/`nregions'{;
			mat taidomega`i'[3,4]=((1+delta_mid_reg_reg)*w_nest[`i',3]
			+beta_reg*(1+delta_top)*w[`i',3]
			+g4p3/sbar`i'[4,1]
			+(b4/sbar`i'[4,1])*(delta_mid_reg_reg*w_nest[`i',3]+w[`i',3]*(1+delta_top)*beta_reg))*x_reg[`i',1]*sbar`i'[4,1];
};



forval i=1/`nregions'{;
			mat taidomega`i'[9,8]=((1+delta_mid_light_light)*w_nest[`i',9]
			+beta_light*(1+delta_top)*w[`i',9]
			+g23p4/sbar`i'[8,1]
			+(b23/sbar`i'[8,1])*(delta_mid_light_light*w_nest[`i',9]+w[`i',9]*(1+delta_top)*beta_light))*x_light[`i',1]*sbar`i'[8,1];
};



forval i=1/`nregions'{;
	mat taidomega`i'[8,9]=((1+delta_mid_light_light)*w_nest[`i',8]
			+beta_light*(1+delta_top)*w[`i',8]
			+g24p3/sbar`i'[9,1]
			+(b24/sbar`i'[9,1])*(delta_mid_light_light*w_nest[`i',8]+w[`i',8]*(1+delta_top)*beta_light))*x_light[`i',1]*sbar`i'[9,1];
};





#delimit;
/*solves for post-merger mark-ups*/
forval i=1/`nregions'{;
    mat postaidomega`i'=aidomega`i'+taidomega`i';
};
forval i=1/`nregions'{;
    mat ivpostaidomega`i'=ivaidomega`i'+ivtaidomega`i';
};

forval i=1/`nregions'{;
    mat paidmu`i'=-inv(postaidomega`i')*X`i'*sbar`i';
};

forval i=1/`nregions'{;
    mat ivpaidmu`i'=-inv(ivpostaidomega`i')*X`i'*sbar`i';
};

/*puts post-merger mark-ups into a matrix with rows corresponding to markets and columns to products*/

mat paidmu=[paidmu1];
mat ivpaidmu=[ivpaidmu1];

forval x=2/`nregions'{;
    mat paidmu=[paidmu,paidmu`x'];
};

mat paidmu=paidmu';


forval x=2/`nregions'{;
    mat ivpaidmu=[ivpaidmu,ivpaidmu`x'];
};

mat ivpaidmu=ivpaidmu';


# delimit;



/*now go back to mata to solve for post-merger prices.  will then calculate percentage increases over 
pre-merger average prices for each region, and then take the average over regions for to yield percentage
price increases for each product.  these average prices are the mxprice variables and the average costs are the
mxcost variables*/

mata;
    /*tsale=st_matrix("tsale");*/
    paidmu=st_matrix("paidmu");
    ivpaidmu=st_matrix("ivpaidmu");
    invaidp=(J(`nregions',`nbrands',1)-paidmu):/aidcost;
    ivinvaidp=(J(`nregions',`nbrands',1)-ivpaidmu):/ivaidcost;
    aidp=((J(`nregions',`nbrands',1):/invaidp)-pbar):/pbar;
    ivaidp=((J(`nregions',`nbrands',1):/ivinvaidp)-pbar):/pbar;
    meanaidcost=mean(aidcost,1);
    ivmeanaidcost=mean(ivaidcost,1);
    meanaidprice=mean(aidp,1);
    ivmeanaidprice=mean(ivaidp,1);
    medianaidprice=mm_median(aidp);
    ivmedianaidprice=mm_median(ivaidp);
end;    



/****************SOLVES UNDER LINEAR DEMAND*****************************/


#delimit cr
mata:
    B_reg=st_matrix("B_reg")
	B_light=st_matrix("B_light")
    V=(st_matrix("vbar1"),st_matrix("vbar2"),st_matrix("vbar3"),st_matrix("vbar4"),st_matrix("vbar5"),st_matrix("vbar6"),st_matrix("vbar7"),st_matrix("vbar8"),st_matrix("vbar9"),st_matrix("vbar10"),
	st_matrix("vbar11"),st_matrix("vbar12"),st_matrix("vbar13"),st_matrix("vbar14"),st_matrix("vbar15"),st_matrix("vbar16"),st_matrix("vbar17"),st_matrix("vbar18"),st_matrix("vbar19"),st_matrix("vbar20"),
	st_matrix("vbar21"),st_matrix("vbar22"),st_matrix("vbar23"),st_matrix("vbar24"),st_matrix("vbar25"),st_matrix("vbar26"),st_matrix("vbar27"),st_matrix("vbar28"),st_matrix("vbar29"),st_matrix("vbar30"),
	st_matrix("vbar31"),st_matrix("vbar32"),st_matrix("vbar33"),st_matrix("vbar34"),st_matrix("vbar35"),st_matrix("vbar36"),st_matrix("vbar37"),st_matrix("vbar38"),st_matrix("vbar39"),st_matrix("vbar40"),
	st_matrix("vbar41"),st_matrix("vbar42"),st_matrix("vbar43"),st_matrix("vbar44"),st_matrix("vbar45"),st_matrix("vbar46"),st_matrix("vbar47"))
	theta_reg=st_matrix("theta_reg")
	theta_light=st_matrix("theta_light")
	ldelta_mid=st_matrix("ldeltamid")
	lbeta_mid=st_matrix("lbetamid")
	ldelta_top=st_matrix("ldelta_top")
	lbeta_reg=st_matrix("lbeta_reg")
		
	
	a=J(9,47,0) 
    c=J(9,47,0)
    pbar=pbar'
	/*create intercepts for regular and then light*/
    for(i=1;i<=47;i++){
        a[1..5,i]=V[1..5,i]-B_reg*pbar[1..5,i]       
    }
	for(i=1;i<=47;i++){
        a[6..9,i]=V[6..9,i]-B_light*pbar[6..9,i]       
    }
	
	B=J(9,9,0)
	B[1..5,1..5]=B_reg
	B[6..9,6..9]=B_light
	/*Create slope matrix here, accounting for middle and upper levels of demand*/
	
	
	mc=J(47,9,0)
	lfullp=J(47,9,0)
	for(j=1;j<=47;j++){
		D_pre=J(9,9,0)
		D_post=J(9,9,0)
		temp=J(9,9,0)
		for(i=1;i<=5;i++){
			D_pre[i,i]=B_reg[i,i]+theta_reg[i,1]*(ldelta_mid[1,1]*w_nest[j,i]+w[j,i]*ldelta_top*lbeta_mid[1,1])
		}
		for(i=6;i<=9;i++){
			D_pre[i,i]=B_light[i-5,i-5]+theta_light[i-5,1]*(ldelta_mid[2,2]*w_nest[j,i]+w[j,i]*ldelta_top*lbeta_mid[2,1])
		}
		D_pre[1,6]=theta_light[1,1]*(ldelta_mid[2,1]*w_nest[j,1]+w[j,1]*ldelta_top*lbeta_mid[2,1])
		D_pre[6,1]=theta_reg[1,1]*(ldelta_mid[1,2]*w_nest[j,6]+w[j,6]*ldelta_top*lbeta_mid[1,1])
		D_pre[2,7]=theta_light[2,1]*(ldelta_mid[2,1]*w_nest[j,2]+w[j,2]*ldelta_top*lbeta_mid[2,1])
		D_pre[7,2]=theta_reg[2,1]*(ldelta_mid[1,2]*w_nest[j,7]+w[j,7]*ldelta_top*lbeta_mid[1,1])
		D_pre[3,8]=theta_light[3,1]*(ldelta_mid[2,1]*w_nest[j,3]+w[j,3]*ldelta_top*lbeta_mid[2,1])
		D_pre[8,3]=theta_reg[3,1]*(ldelta_mid[1,2]*w_nest[j,8]+w[j,8]*ldelta_top*lbeta_mid[1,1])
		D_pre[4,9]=theta_light[4,1]*(ldelta_mid[2,1]*w_nest[j,4]+w[j,4]*ldelta_top*lbeta_mid[2,1])
		D_pre[9,4]=theta_reg[4,1]*(ldelta_mid[1,2]*w_nest[j,9]+w[j,9]*ldelta_top*lbeta_mid[1,1])
		
		temp[3,4]=B_reg[4,3]+theta_reg[4,1]*(ldelta_mid[1,1]*w_nest[j,3]+w[j,3]*ldelta_top*lbeta_mid[1,1])
		temp[4,3]=B_reg[3,4]+theta_reg[3,1]*(ldelta_mid[1,1]*w_nest[j,4]+w[j,4]*ldelta_top*lbeta_mid[1,1])
		temp[8,9]=B_light[4,3]+theta_light[4,1]*(ldelta_mid[2,2]*w_nest[j,8]+w[j,8]*ldelta_top*lbeta_mid[2,1])
		temp[9,8]=B_light[3,4]+theta_light[3,1]*(ldelta_mid[2,2]*w_nest[j,9]+w[j,9]*ldelta_top*lbeta_mid[2,1])
		
		temp[3,9]=theta_light[3,1]*(ldelta_mid[2,1]*w_nest[j,3]+w[j,3]*ldelta_top*lbeta_mid[2,1])
		temp[9,3]=theta_reg[3,1]*(ldelta_mid[1,2]*w_nest[j,9]+w[j,9]*ldelta_top*lbeta_mid[1,1])
		temp[4,8]=theta_light[4,1]*(ldelta_mid[2,1]*w_nest[j,4]+w[j,4]*ldelta_top*lbeta_mid[2,1])
		temp[8,4]=theta_reg[4,1]*(ldelta_mid[1,2]*w_nest[j,8]+w[j,8]*ldelta_top*lbeta_mid[1,1])
		
		D_post=temp+D_pre
		
		mcind=luinv(D_pre)*(a[.,j]+B*pbar[.,j]+D_pre*pbar[.,j])
		mc[j,.]=mcind'
		pind=luinv(B+D_post)*(D_post*mc[j,.]'-a[.,j])
		lfullp[j,.]=pind'
		linans=mm_median((lfullp-pbar'):/pbar')
	}	
	
end

mata:
    ivB_reg=st_matrix("ivB_reg")
	ivB_light=st_matrix("ivB_light")
    V=(st_matrix("vbar1"),st_matrix("vbar2"),st_matrix("vbar3"),st_matrix("vbar4"),st_matrix("vbar5"),st_matrix("vbar6"),st_matrix("vbar7"),st_matrix("vbar8"),st_matrix("vbar9"),st_matrix("vbar10"),
	st_matrix("vbar11"),st_matrix("vbar12"),st_matrix("vbar13"),st_matrix("vbar14"),st_matrix("vbar15"),st_matrix("vbar16"),st_matrix("vbar17"),st_matrix("vbar18"),st_matrix("vbar19"),st_matrix("vbar20"),
	st_matrix("vbar21"),st_matrix("vbar22"),st_matrix("vbar23"),st_matrix("vbar24"),st_matrix("vbar25"),st_matrix("vbar26"),st_matrix("vbar27"),st_matrix("vbar28"),st_matrix("vbar29"),st_matrix("vbar30"),
	st_matrix("vbar31"),st_matrix("vbar32"),st_matrix("vbar33"),st_matrix("vbar34"),st_matrix("vbar35"),st_matrix("vbar36"),st_matrix("vbar37"),st_matrix("vbar38"),st_matrix("vbar39"),st_matrix("vbar40"),
	st_matrix("vbar41"),st_matrix("vbar42"),st_matrix("vbar43"),st_matrix("vbar44"),st_matrix("vbar45"),st_matrix("vbar46"),st_matrix("vbar47"))
	ivtheta_reg=st_matrix("ivtheta_reg")
	ivtheta_light=st_matrix("ivtheta_light")
	ivldelta_mid=st_matrix("livdeltamid")
	ivlbeta_mid=st_matrix("livbetamid")
	ivldelta_top=st_matrix("ivldelta_top")
	ivlbeta_reg=st_matrix("livbeta_reg")
		
	
	a=J(9,47,0) 
    c=J(9,47,0)
    
	/*create intercepts for regular and then light*/
    for(i=1;i<=47;i++){
        a[1..5,i]=V[1..5,i]-ivB_reg*pbar[1..5,i]       
    }
	for(i=1;i<=47;i++){
        a[6..9,i]=V[6..9,i]-ivB_light*pbar[6..9,i]       
    }
	
	ivB=J(9,9,0)
	ivB[1..5,1..5]=ivB_reg
	ivB[6..9,6..9]=ivB_light
	/*Create slope matrix here, accounting for middle and upper levels of demand*/
	
	
	ivmc=J(47,9,0)
	ivlfullp=J(47,9,0)
	for(j=1;j<=47;j++){
		ivD_pre=J(9,9,0)
		ivD_post=J(9,9,0)
		temp=J(9,9,0)
		for(i=1;i<=5;i++){
			ivD_pre[i,i]=ivB_reg[i,i]+ivtheta_reg[i,1]*(ivldelta_mid[1,1]*w_nest[j,i]+w[j,i]*ivldelta_top*ivlbeta_mid[1,1])
		}
		for(i=6;i<=9;i++){
			ivD_pre[i,i]=ivB_light[i-5,i-5]+ivtheta_light[i-5,1]*(ivldelta_mid[2,2]*w_nest[j,i]+w[j,i]*ivldelta_top*ivlbeta_mid[2,1])
		}
		ivD_pre[1,6]=ivtheta_light[1,1]*(ivldelta_mid[2,1]*w_nest[j,1]+w[j,1]*ivldelta_top*ivlbeta_mid[2,1])
		ivD_pre[6,1]=ivtheta_reg[1,1]*(ivldelta_mid[1, 2]*w_nest[j,6]+w[j,6]*ivldelta_top*ivlbeta_mid[1,1])
		ivD_pre[2,7]=ivtheta_light[2,1]*(ivldelta_mid[2,1]*w_nest[j,2]+w[j,2]*ivldelta_top*ivlbeta_mid[2,1])
		ivD_pre[7,2]=ivtheta_reg[2,1]*(ivldelta_mid[1,2]*w_nest[j,7]+w[j,7]*ivldelta_top*ivlbeta_mid[1,1])
		ivD_pre[3,8]=ivtheta_light[3,1]*(ivldelta_mid[2,1]*w_nest[j,3]+w[j,3]*ivldelta_top*ivlbeta_mid[2,1])
		ivD_pre[8,3]=ivtheta_reg[3,1]*(ivldelta_mid[1,2]*w_nest[j,8]+w[j,8]*ivldelta_top*ivlbeta_mid[1,1])
		ivD_pre[4,9]=ivtheta_light[4,1]*(ivldelta_mid[2,1]*w_nest[j,4]+w[j,4]*ivldelta_top*ivlbeta_mid[2,1])
		ivD_pre[9,4]=ivtheta_reg[4,1]*(ivldelta_mid[1,2]*w_nest[j,9]+w[j,9]*ivldelta_top*ivlbeta_mid[1,1])
		
		temp[3,4]=ivB_reg[4,3]+ivtheta_reg[4,1]*(ivldelta_mid[1,1]*w_nest[j,3]+w[j,3]*ivldelta_top*ivlbeta_mid[1,1])
		temp[4,3]=ivB_reg[3,4]+ivtheta_reg[3,1]*(ivldelta_mid[1,1]*w_nest[j,4]+w[j,4]*ivldelta_top*ivlbeta_mid[1,1])
		temp[8,9]=ivB_light[4,3]+ivtheta_light[4,1]*(ivldelta_mid[2,2]*w_nest[j,8]+w[j,8]*ivldelta_top*ivlbeta_mid[2,1])
		temp[9,8]=ivB_light[3,4]+ivtheta_light[3,1]*(ivldelta_mid[2,2]*w_nest[j,9]+w[j,9]*ivldelta_top*ivlbeta_mid[2,1])
		
		temp[3,9]=ivtheta_light[3,1]*(ivldelta_mid[2,1]*w_nest[j,3]+w[j,3]*ivldelta_top*ivlbeta_mid[2,1])
		temp[9,3]=ivtheta_reg[3,1]*(ivldelta_mid[1,2]*w_nest[j,9]+w[j,9]*ivldelta_top*ivlbeta_mid[1,1])
		temp[4,8]=ivtheta_light[4,1]*(ivldelta_mid[2,1]*w_nest[j,4]+w[j,4]*ivldelta_top*ivlbeta_mid[2,1])
		temp[8,4]=ivtheta_reg[4,1]*(ivldelta_mid[1,2]*w_nest[j,8]+w[j,8]*ivldelta_top*ivlbeta_mid[1,1])
		
		ivD_post=temp+ivD_pre
		
		ivmcind=luinv(ivD_pre)*(a[.,j]+ivB*pbar[.,j]+ivD_pre*pbar[.,j])
		ivmc[j,.]=ivmcind'
		ivpind=luinv(ivB+ivD_post)*(ivD_post*ivmc[j,.]'-a[.,j])
		ivlfullp[j,.]=ivpind'
		ivlinans=mm_median((ivlfullp-pbar'):/pbar')
	}	
end

/*this solves for the exact equilibrium for each model*/
#delimit cr
local nbrands=9
local nregions=47
mata:
	pbar=pbar'
	aidfullp=J(`nregions',`nbrands',0)
	check=J(`nregions',1,0)
    for(i=1;i<=`nregions';i++){
		a=csolve(&msyraids(),pbar[i,.]',1e-6,400,sbar[i,.]',pbar[i,.]',aidcost[i,.]',g_reg,g_light,b_reg,b_light,w[i,.]',w_nest[i,.]',X[i,.]',delta,beta,delta_top)		
        aidfullp[i,.]=a[1..`nbrands',.]'
        ws=msyraids(a[1..9,.],sbar[i,.]',pbar[i,.]',aidcost[i,.]',g_reg,g_light,b_reg,b_light,w[i,.]',w_nest[i,.]',X[i,.]',delta,beta,delta_top)		
		check[i,1]=((a[10,1]+(a[1..rows(a)-1,.]==J(rows(a)-1,1,.))+(ws[1,.]==.)+(ws[2,.]==.)+(ws[3,.]==.)+(ws[4,.]==.)+(ws[5,.]==.)+(ws[6,.]==.)+(ws[7,.]==.)+(ws[8,.]==.)+(ws[9,.]==.))==0)
    }
    ans=(aidfullp-pbar):/pbar
    faids=mm_median(select(ans,check))    	
end


mata:
	ivaidfullp=J(`nregions',`nbrands',0)
	ivcheck=J(`nregions',1,0)
    for(i=1;i<=`nregions';i++){
		a=csolve(&msyraids(),pbar[i,.]',1e-6,400,sbar[i,.]',pbar[i,.]',ivaidcost[i,.]',ivg_reg,ivg_light,ivb_reg,ivb_light,w[i,.]',w_nest[i,.]',X[i,.]',ivdelta,ivbeta,ivdelta_top)		
        ivaidfullp[i,.]=a[1..`nbrands',.]'
		ws=msyraids(a[1..9,.],sbar[i,.]',pbar[i,.]',ivaidcost[i,.]',ivg_reg,ivg_light,ivb_reg,ivb_light,w[i,.]',w_nest[i,.]',X[i,.]',ivdelta,ivbeta,ivdelta_top)
		ivcheck[i,1]=((a[10,1]+(a[1..rows(a)-1,.]==J(rows(a)-1,1,.))+(ws[1,.]==.)+(ws[2,.]==.)+(ws[3,.]==.)+(ws[4,.]==.)+(ws[5,.]==.)+(ws[6,.]==.)+(ws[7,.]==.)+(ws[8,.]==.)+(ws[9,.]==.)  )==0)		
    }	
    ivans=(ivaidfullp-pbar):/pbar
    ivfaids=mm_median(select(ivans,ivcheck))    	
end

mata:
	linans
	ivlinans
	faids
	ivfaids
end:	




/*
do msyrcosts
mata:
    mata matuse post
    cpbar=pbar
    cpbar[.,3]=pbarpost[.,3]		
    cpbar[.,4]=pbarpost[.,4]

    caidfullp=J(`nbrands',`nregions',0)
    for(i=1;i<=cols(pbar');i++){
        pinit=cpbar[i,.]'
        ca=csolve(&msyrcosts(),pinit,1e-6,400,sbar[i,.]',cpbar[i,.]',aidcost[i,.]',g,b,w[i,.]',tsale[i])
        ca=ca[1..`nbrands',.]
        caidfullp[.,i]=ca
    }
	caidfullp=caidfullp'
	mean((caidfullp[.,3]-aidcost[.,3]):/aidcost[.,3])
	mean((caidfullp[.,4]-aidcost[.,4]):/aidcost[.,4])    	
end
*/	


#delimit;
save postsyrstart, replace; 
drop _all;


#delimit;
/***Now that we have our point estimates of costs and approximate simulated price changes, we draw from asymptotic approx of distributions
of estimators and recalculate for each draw.  will look at empirical distribution of these calculated costs and changes to do inference.****/
drop _all;


/*number of bootstrap iterations*/
local its=1000;
local aids_reg "laidpreg1 laidpreg2 laidpreg3 laidpreg4 `reglist1' `Mlist1' `regionxbrand1' bra11 bra12 bra13 bra14 count11 count12 count13 count14";
local aids_light "laidplight1 laidplight2 laidplight3 `reglist2' `Mlist2' `regionxbrand2' bra21 bra22 bra23 count21 count22 count23";
local aids_mid "lexp1 lexp2 lp1_1 lp1_2 lp2_1 lp2_2 `femid' `Mlistmid' nest1 nest2 count1 count2";
local aids_top "indexall month2 month3 month4 month5 month6 month7 month8 month9 month10 month11 month12 count reg2 reg3 reg4 reg5 reg6 reg7 reg8 reg9 reg10 reg11 reg12 reg13 reg14 reg15 reg16 reg17 reg18 reg19 reg20 reg21 reg22 reg23 reg24 reg25 reg26 reg27 reg28 reg29 reg30 reg31 reg32 reg33 reg34 reg35 reg36 reg37 reg38 reg39 reg40 reg41 reg42 reg43 reg44 reg45 reg46 reg47 const";

 

drawnorm `aids_reg', n(`its') mean(ivaidB) cov(ivaidVC);
drawnorm `aids_light', n(`its') mean(ivaidB2) cov(ivaidVC2); 
drawnorm `aids_mid', n(`its') mean(ivaidMidB) cov(ivaidMidVC);
drawnorm `aids_top', n(`its') mean(iv_b_top) cov(iv_vc_top);


#delimit;
foreach x of local list1{;
    gen lpreg_5_`x'=-lpreg_1_`x'-lpreg_2_`x'-lpreg_3_`x'-lpreg_4_`x';
};
gen laidpreg5=-laidpreg1-laidpreg2-laidpreg3-laidpreg4;

forval x=1/5{;
	rename lpreg_`x'_AUNT_JEMIMA lpreg_`x'_1;
	rename lpreg_`x'_HUNGRY_JACK lpreg_`x'_2;
	rename lpreg_`x'_LOG_CABIN lpreg_`x'_3;
	rename lpreg_`x'_MRS_BUTTERWORTH lpreg_`x'_4;
	rename lpreg_`x'_PRIVATE_LABEL lpreg_`x'_5;
};

foreach x of local list2{;
    gen lplight_4_`x'=-lplight_1_`x'-lplight_2_`x'-lplight_3_`x';
};
gen laidplight4=-laidplight1-laidplight2-laidplight3;


forval x=1/4{;
	rename lplight_`x'_AUNT_JEMIMA_LIGHT lplight_`x'_1;
	rename lplight_`x'_JACK_LITE lplight_`x'_2;
	rename lplight_`x'_LOG_CABIN_LIGHT lplight_`x'_3;
	rename lplight_`x'_MRS_B_LITE lplight_`x'_4;
};

drop M11b1-count14 M21b1-count23 fe1_2-nest2 count1 count2  month2-reg47;

mkmat _all;


use postsyrstart, replace;
keep if merged==1;


/*bsbarx is a 7x`its' matrix with columns equal to the means of revenue shares taken over bootstrap subsamples
of region x.  bpbarx, bvbbarx defined similarly.*/

forval x=1/`nregions'{;
    mat bsbar`x'=J(9,`its',0);
    mat bpbar`x'=J(9,`its',0);
	mat x_reg`x'=J(1,`its',0);
	mat x_light`x'=J(1,`its',0);
};

forval x=1/`nregions'{;
    forval w=1/`its'{;
        preserve;
            keep if reg`x'==1;
            bsample;
			qui sum expmid if reg`x'==1&nest==1;
			mat x_reg`x'[1,`w']=_result(3);
			qui sum expmid if reg`x'==1&nest==2;
			mat x_light`x'[1,`w']=_result(3);
			
			forval y=1/`nbrands'{;
                qui sum share if bra`y'==1;
                matrix bsbar`x'[`y',`w']=_result(3);
                qui sum price if bra`y'==1;
                matrix bpbar`x'[`y',`w']=_result(3);
            };
        restore;
    };
};

#delimit;
local its=1000;
local nbrands=9;
local nregions=47;
local nbrands1=5;
local nbrands2=4;
mata;
    bsaidcost=J(`its',9,0);
    bsaidp=J(`its',9,0);
end;

#delimit;
local its=1000;
forval z=1/`its'{;
    /*makes pre-merger "ownership" matrices*/
	forval i=1/`nregions'{;
        mat define aidomega`i'=J(`nbrands',`nbrands',0);
    };

	
	
	forval i=1/`nregions'{;
		mat define X`i'=J(`nbrands',`nbrands',0);
	};

	
	forval i=1/`nregions'{;
		forval x=1/`nbrands1'{;
			mat aidomega`i'[`x',`x']=(-1+(1+lp1_1[`z',1])*w_nest[`i',`x']
			+lexp1[`z',1]*(1+indexall[`z',1])*w[`i',`x']
			+lpreg_`x'_`x'[`z',1]/bsbar`i'[`x',`z']
			+(laidpreg`x'[`z',1]/bsbar`i'[`x',`z'])*(lp1_1[`z',1]*w_nest[`i',`x']+w[`i',`x']*(1+indexall[`z',1])*lexp1[`z',1]))*x_reg`i'[1,`z']*bsbar`i'[`x',`z'];
		};
	};

	
	
	forval i=1/`nregions'{;
		forval x=1/`nbrands2'{;
			mat aidomega`i'[5+`x',5+`x']=(-1+(1+lp2_2[`z',1])*w_nest[`i',5+`x']
			+lexp2[`z',1]*(1+indexall[`z',1])*w[`i',5+`x']
			+lplight_`x'_`x'[`z',1]/bsbar`i'[5+`x',`z']
			+(laidplight`x'[`z',1]/bsbar`i'[5+`x',`z'])*(lp2_2[`z',1]*w_nest[`i',5+`x']+w[`i',5+`x']*(1+indexall[`z',1])*lexp2[`z',1]))*x_light`i'[1,`z']*bsbar`i'[`x'+5,`z'];
		};
	};

	
	forval i=1/`nregions'{;
		mat aidomega`i'[6,1]=((1+laidpreg1[`z',1]/bsbar`i'[1,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',6]+lp1_2[`z',1]*w_nest[`i',6]))*x_reg`i'[1,`z']*bsbar`i'[1,`z'];
	};
	
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[1,6]=((1+laidplight1[`z',1]/bsbar`i'[6,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',1]+lp2_1[`z',1]*w_nest[`i',1]))*x_light`i'[1,`z']*bsbar`i'[6,`z'];
	};
	
	
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[7,2]=((1+laidpreg2[`z',1]/bsbar`i'[2,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',7]+lp1_2[`z',1]*w_nest[`i',7]))*x_reg`i'[1,`z']*bsbar`i'[2,`z'];
	};
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[2,7]=((1+laidplight2[`z',1]/bsbar`i'[7,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',2]+lp2_1[`z',1]*w_nest[`i',2]))*x_light`i'[1,`z']*bsbar`i'[7,`z'];
	};
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[8,3]=((1+laidpreg3[`z',1]/bsbar`i'[3,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',8]+lp1_2[`z',1]*w_nest[`i',8]))*x_reg`i'[1,`z']*bsbar`i'[3,`z'];
	};
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[3,8]=((1+laidplight3[`z',1]/bsbar`i'[8,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',3]+lp2_1[`z',1]*w_nest[`i',3]))*x_light`i'[1,`z']*bsbar`i'[8,`z'];
	};
	
	forval i=1/`nregions'{;
		mat aidomega`i'[9,4]=((1+laidpreg4[`z',1]/bsbar`i'[4,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',9]+lp1_2[`z',1]*w_nest[`i',9]))*x_reg`i'[1,`z']*bsbar`i'[4,`z'];
	};
	
	
	forval i=1/`nregions'{;
		mat aidomega`i'[4,9]=((1+laidplight4[`z',1]/bsbar`i'[9,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',4]+lp2_1[`z',1]*w_nest[`i',4]))*x_light`i'[1,`z']*bsbar`i'[9,`z'];
	};
	
	
	forval i=1/`nregions'{;
		forval x=1/`nbrands1'{;
			mat X`i'[`x',`x']=x_reg`i'[1,`z'];
		};
	};

	forval i=1/`nregions'{;
		forval x=1/`nbrands2'{;
			mat X`i'[`x'+5,`x'+5]=x_light`i'[1,`z'];
		};
	};

	forval i=1/`nregions'{;
		mat aidmu`i'=J(`nbrands',1,0);
		mat aidmu`i'=-inv(aidomega`i')*X`i'*bsbar`i'[1..`nbrands',`z'];		
	};



/*mark-up matrix has 10 columns, which correspond to regional markets, and 7 rows corresponding to brands.  
aidmu are p-mc's with aid system, lmu are p=mc's with linear system, llmu are lerner indices with log-linear system
this just puts them in matrix form as described above*/
	mat aidmu=[aidmu1];
	mat pbar=[bpbar1[1..`nbrands',`z']];
    mat sbar=[bsbar1[1..`nbrands',`z']];
	forval x=2/`nregions'{;
		mat aidmu=[aidmu,aidmu`x'];
		mat pbar=[pbar, bpbar`x'[1..`nbrands',`z']];
		mat sbar=[sbar, bsbar`x'[1..`nbrands',`z']];
	};

	mat pbar=pbar';
	/*mat lpbar=lpbar';*/
	mat sbar=sbar';
	mat aidmu=aidmu';

/*get costs by solving mark-ups.  because element-by-element operations are a pain in stata, i do this in mata*/
/*input markup matrices, average prices and shares, and post-merger ownership matrices into stata*/
	mata: aidmu=st_matrix("aidmu");
	mata: pbar=st_matrix("pbar"); 
	mata: sbar=st_matrix("sbar");
	mata: aidcost=pbar:*(J(`nregions',`nbrands',1)-aidmu);
	mata: st_matrix("aidcost", aidcost);
	

	
    /*defines post-merger ownership matrix of elasticities*/
    mat taidomega=J(`nbrands',`nbrands',0);
	forval i=1/`nregions'{;
		mat taidomega`i'[8,4]=((1+laidpreg4[`z',1]/bsbar`i'[4,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',8]+lp1_2[`z',1]*w_nest[`i',8]))*x_reg`i'[1,`z']*bsbar`i'[4,`z'];
	};

	forval i=1/`nregions'{;
		mat taidomega`i'[4,8]=((1+laidplight3[`z',1]/bsbar`i'[8,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',4]+lp2_1[`z',1]*w_nest[`i',4]))*x_light`i'[1,`z']*bsbar`i'[8,`z'];
	};

	forval i=1/`nregions'{;
		mat taidomega`i'[9,3]=((1+laidpreg3[`z',1]/bsbar`i'[3,`z'])*(lexp1[`z',1]*(1+indexall[`z',1])*w[`i',9]+lp1_2[`z',1]*w_nest[`i',9]))*x_reg`i'[1,`z']*bsbar`i'[3,`z'];
	};


	forval i=1/`nregions'{;
		mat taidomega`i'[3,9]=((1+laidplight4[`z',1]/bsbar`i'[9,`z'])*(lexp2[`z',1]*(1+indexall[`z',1])*w[`i',3]+lp2_1[`z',1]*w_nest[`i',3]))*x_light`i'[1,`z']*bsbar`i'[9,`z'];
	};



	forval i=1/`nregions'{;            
			mat taidomega`i'[4,3]=((1+lp1_1[`z',1])*w_nest[`i',4]
			+lexp1[`z',1]*(1+indexall[`z',1])*w[`i',4]
			+lpreg_4_3[`z',1]/bsbar`i'[3,`z']
			+(laidpreg3[`z',1]/bsbar`i'[3,`z'])*(lp1_1[`z',1]*w_nest[`i',4]+w[`i',4]*(1+indexall[`z',1])*lexp1[`z',1]))*x_reg`i'[1,`z']*bsbar`i'[3,`z'];
	};

	
	forval i=1/`nregions'{;
			mat taidomega`i'[3,4]=((1+lp1_1[`z',1])*w_nest[`i',3]
			+lexp1[`z',1]*(1+indexall[`z',1])*w[`i',3]
			+lpreg_3_4[`z',1]/bsbar`i'[4,`z']
			+(laidpreg4[`z',1]/bsbar`i'[4,`z'])*(lp1_1[`z',1]*w_nest[`i',3]+w[`i',3]*(1+indexall[`z',1])*lexp1[`z',1]))*x_reg`i'[1,`z']*bsbar`i'[4,`z'];
	};



	forval i=1/`nregions'{;
			mat taidomega`i'[9,8]=((1+lp2_2[`z',1])*w_nest[`i',9]
			+lexp2[`z',1]*(1+indexall[`z',1])*w[`i',9]
			+lplight_4_3[`z',1]/bsbar`i'[8,`z']
			+(laidplight3[`z',1]/bsbar`i'[8,`z'])*(lp2_2[`z',1]*w_nest[`i',9]+w[`i',9]*(1+indexall[`z',1])*lexp2[`z',1]))*x_light`i'[1,`z']*bsbar`i'[8,`z'];
	};



	forval i=1/`nregions'{;
		mat taidomega`i'[8,9]=((1+lp2_2[`z',1])*w_nest[`i',8]
			+lexp2[`z',1]*(1+indexall[`z',1])*w[`i',8]
			+lplight_3_4[`z',1]/bsbar`i'[9,`z']
			+(laidplight4[`z',1]/bsbar`i'[9,`z'])*(lp2_2[`z',1]*w_nest[`i',8]+w[`i',8]*(1+indexall[`z',1])*lexp2[`z',1]))*x_light`i'[1,`z']*bsbar`i'[9,`z'];
	};





	#delimit;
	/*solves for post-merger mark-ups*/
	forval i=1/`nregions'{;
		mat postaidomega`i'=aidomega`i'+taidomega`i';
	};


    

	/*solves for post-merger mark-ups*/
    forval i=1/`nregions'{;
        mat define paidmu`i'=J(`nbrands',1,0);
        mat paidmu`i'=-inv(postaidomega`i')*X`i'*bsbar`i'[1..`nbrands',`z'];
    };
    /*puts post-merger mark-ups into a matrix with rows corresponding to markets and columns to products*/
    mat paidmu=[paidmu1];
    forval x=2/`nregions'{;
        mat paidmu=[paidmu,paidmu`x'];
    };
    mat paidmu=paidmu';

    /*now go back to mata to solve for post-merger prices.  will then calculate percentage increases over
    pre-merger average prices for each region, and then take the average over regions for to yield percentage
    price increases for each product.  these average prices are the mxprice variables and the average costs are the
    mxcost variables*/

    mata: paidmu=st_matrix("paidmu");
    mata: invaidp=(J(47,9,1)-paidmu):/aidcost;
    mata: aidp=((J(47,9,1):/invaidp)-pbar):/pbar;
    mata: maidcost=mean(aidcost,1);
    mata: maidprice=mean(aidp,1);
    mata: bsaidcost[`z',.]=maidcost;
    mata: bsaidp[`z',.]=maidprice;
    mata: st_matrix("bsaidcost",bsaidcost);
    mata: st_matrix("bsaidp", bsaidp);
};






/*save nowait*/
#delimit;
drop _all;
svmat bsaidcost;
svmat bsaidp;
save sampledists, replace;

su bsaidcost*, detail;
su bsaidp3, detail;
su bsaidp4, detail;
su bsaidp8, detail;
su bsaidp9, detail;

#delimit;
/*twoway hist bsaidp3 if bsaidp4<1 & bsaidp4>0, title("AID") subtitle("Pennzoil") xtitle("% Price Change in Hundredths") saving(pennapprox, replace);
twoway hist bsaidp4 if bsaidp4<1 & bsaidp4>0, title("AID") subtitle("QuakerState") xtitle("% Price Change in Hundredths") saving(quakerapprox, replace);
*/

#delimit cr;
set more off
mata:
    bp1=st_matrix("bpbar1")
    bp2=st_matrix("bpbar2")
    bp3=st_matrix("bpbar3")
    bp4=st_matrix("bpbar4")
    bp5=st_matrix("bpbar5")
    bp6=st_matrix("bpbar6")
    bp7=st_matrix("bpbar7")
    bp8=st_matrix("bpbar8")
    bp9=st_matrix("bpbar9")
    bp10=st_matrix("bpbar10")
    bp11=st_matrix("bpbar11")
    bp12=st_matrix("bpbar12")
    bp13=st_matrix("bpbar13")
    bp14=st_matrix("bpbar14")
    bp15=st_matrix("bpbar15")
    bp16=st_matrix("bpbar16")
    bp17=st_matrix("bpbar17")
    bp18=st_matrix("bpbar18")
    bp19=st_matrix("bpbar19")
    bp20=st_matrix("bpbar20")
    bp21=st_matrix("bpbar21")
    bp22=st_matrix("bpbar22")
    bp23=st_matrix("bpbar23")
    bp24=st_matrix("bpbar24")
    bp25=st_matrix("bpbar25")
    bp26=st_matrix("bpbar26")
    bp27=st_matrix("bpbar27")
    bp28=st_matrix("bpbar28")
    bp29=st_matrix("bpbar29")
    bp30=st_matrix("bpbar30")
    bp31=st_matrix("bpbar31")
    bp32=st_matrix("bpbar32")
    bp33=st_matrix("bpbar33")
    bp34=st_matrix("bpbar34")
    bp35=st_matrix("bpbar35")
    bp36=st_matrix("bpbar36")
    bp37=st_matrix("bpbar37")
    bp38=st_matrix("bpbar38")
    bp39=st_matrix("bpbar39")
    bp40=st_matrix("bpbar40")
    bp41=st_matrix("bpbar41")
    bp42=st_matrix("bpbar42")
    bp43=st_matrix("bpbar43")
    bp44=st_matrix("bpbar44")
    bp45=st_matrix("bpbar45")
    bp46=st_matrix("bpbar46")
    bp47=st_matrix("bpbar47")

    ppoint=J(1,47,NULL)
    ppoint[1]=&bp1
    ppoint[2]=&bp2
    ppoint[3]=&bp3
    ppoint[4]=&bp4
    ppoint[5]=&bp5
    ppoint[6]=&bp6
    ppoint[7]=&bp7
    ppoint[8]=&bp8
    ppoint[9]=&bp9
    ppoint[10]=&bp10
    ppoint[11]=&bp11
    ppoint[12]=&bp12
    ppoint[13]=&bp13
    ppoint[14]=&bp14
    ppoint[15]=&bp15
    ppoint[16]=&bp16
    ppoint[17]=&bp17
    ppoint[18]=&bp18
    ppoint[19]=&bp19
    ppoint[20]=&bp20
    ppoint[21]=&bp21
    ppoint[22]=&bp22
    ppoint[23]=&bp23
    ppoint[24]=&bp24
    ppoint[25]=&bp25
    ppoint[26]=&bp26
    ppoint[27]=&bp27
    ppoint[28]=&bp28
    ppoint[29]=&bp29
    ppoint[30]=&bp30
    ppoint[31]=&bp31
    ppoint[32]=&bp32
    ppoint[33]=&bp33
    ppoint[34]=&bp34
    ppoint[35]=&bp35
    ppoint[36]=&bp36
    ppoint[37]=&bp37
    ppoint[38]=&bp38
    ppoint[39]=&bp39
    ppoint[40]=&bp40
    ppoint[41]=&bp41
    ppoint[42]=&bp42
    ppoint[43]=&bp43
    ppoint[44]=&bp44
    ppoint[45]=&bp45
    ppoint[46]=&bp46
    ppoint[47]=&bp47
    

        
    bs1=st_matrix("bsbar1")
    bs2=st_matrix("bsbar2")
    bs3=st_matrix("bsbar3")
    bs4=st_matrix("bsbar4")
    bs5=st_matrix("bsbar5")
    bs6=st_matrix("bsbar6")
    bs7=st_matrix("bsbar7")
    bs8=st_matrix("bsbar8")
    bs9=st_matrix("bsbar9")
    bs10=st_matrix("bsbar10")
    bs11=st_matrix("bsbar11")
    bs12=st_matrix("bsbar12")
    bs13=st_matrix("bsbar13")
    bs14=st_matrix("bsbar14")
    bs15=st_matrix("bsbar15")
    bs16=st_matrix("bsbar16")
    bs17=st_matrix("bsbar17")
    bs18=st_matrix("bsbar18")
    bs19=st_matrix("bsbar19")
    bs20=st_matrix("bsbar20")
    bs21=st_matrix("bsbar21")
    bs22=st_matrix("bsbar22")
    bs23=st_matrix("bsbar23")
    bs24=st_matrix("bsbar24")
    bs25=st_matrix("bsbar25")
    bs26=st_matrix("bsbar26")
    bs27=st_matrix("bsbar27")
    bs28=st_matrix("bsbar28")
    bs29=st_matrix("bsbar29")
    bs30=st_matrix("bsbar30")
    bs31=st_matrix("bsbar31")
    bs32=st_matrix("bsbar32")
    bs33=st_matrix("bsbar33")
    bs34=st_matrix("bsbar34")
    bs35=st_matrix("bsbar35")
    bs36=st_matrix("bsbar36")
    bs37=st_matrix("bsbar37")
    bs38=st_matrix("bsbar38")
    bs39=st_matrix("bsbar39")
    bs40=st_matrix("bsbar40")
    bs41=st_matrix("bsbar41")
    bs42=st_matrix("bsbar42")
    bs43=st_matrix("bsbar43")
    bs44=st_matrix("bsbar44")
    bs45=st_matrix("bsbar45")
    bs46=st_matrix("bsbar46")
    bs47=st_matrix("bsbar47")
    

    spoint=J(1,47,NULL)
    spoint[1]=&bs1
    spoint[2]=&bs2
    spoint[3]=&bs3
    spoint[4]=&bs4
    spoint[5]=&bs5
    spoint[6]=&bs6
    spoint[7]=&bs7
    spoint[8]=&bs8
    spoint[9]=&bs9
    spoint[10]=&bs10
    spoint[11]=&bs11
    spoint[12]=&bs12
    spoint[13]=&bs13
    spoint[14]=&bs14
    spoint[15]=&bs15
    spoint[16]=&bs16
    spoint[17]=&bs17
    spoint[18]=&bs18
    spoint[19]=&bs19
    spoint[20]=&bs20
    spoint[21]=&bs21
    spoint[22]=&bs22
    spoint[23]=&bs23
    spoint[24]=&bs24
    spoint[25]=&bs25
    spoint[26]=&bs26
    spoint[27]=&bs27
    spoint[28]=&bs28
    spoint[29]=&bs29
    spoint[30]=&bs30
    spoint[31]=&bs31
    spoint[32]=&bs32
    spoint[33]=&bs33
    spoint[34]=&bs34
    spoint[35]=&bs35
    spoint[36]=&bs36
    spoint[37]=&bs37
    spoint[38]=&bs38
    spoint[39]=&bs39
    spoint[40]=&bs40
    spoint[41]=&bs41
    spoint[42]=&bs42
    spoint[43]=&bs43
    spoint[44]=&bs44
    spoint[45]=&bs45
    spoint[46]=&bs46
    spoint[47]=&bs47
    

/*
    bv1=st_matrix("bvbar1")
    bv2=st_matrix("bvbar2")
    bv3=st_matrix("bvbar3")
    bv4=st_matrix("bvbar4")
    bv5=st_matrix("bvbar5")
    bv6=st_matrix("bvbar6")
    bv7=st_matrix("bvbar7")
    bv8=st_matrix("bvbar8")
    bv9=st_matrix("bvbar9")
    bv10=st_matrix("bvbar10")
    bv11=st_matrix("bvbar11")
    bv12=st_matrix("bvbar12")
    bv13=st_matrix("bvbar13")
    bv14=st_matrix("bvbar14")
    bv15=st_matrix("bvbar15")
    bv16=st_matrix("bvbar16")
    bv17=st_matrix("bvbar17")
    bv18=st_matrix("bvbar18")
    bv19=st_matrix("bvbar19")
    bv20=st_matrix("bvbar20")
    bv21=st_matrix("bvbar21")
    bv22=st_matrix("bvbar22")
    bv23=st_matrix("bvbar23")
    bv24=st_matrix("bvbar24")
    bv25=st_matrix("bvbar25")
    bv26=st_matrix("bvbar26")
    bv27=st_matrix("bvbar27")
    bv28=st_matrix("bvbar28")
    bv29=st_matrix("bvbar29")
    bv30=st_matrix("bvbar30")
    bv31=st_matrix("bvbar31")
    bv32=st_matrix("bvbar32")
    bv33=st_matrix("bvbar33")
    bv34=st_matrix("bvbar34")
    bv35=st_matrix("bvbar35")
    bv36=st_matrix("bvbar36")
    bv37=st_matrix("bvbar37")
    bv38=st_matrix("bvbar38")
    bv39=st_matrix("bvbar39")
    bv40=st_matrix("bvbar40")

    vpoint=J(1,40,NULL)
    vpoint[1]=&bv1
    vpoint[2]=&bv2
    vpoint[3]=&bv3
    vpoint[4]=&bv4
    vpoint[5]=&bv5
    vpoint[6]=&bv6
    vpoint[7]=&bv7
    vpoint[8]=&bv8
    vpoint[9]=&bv9
    vpoint[10]=&bv10
    vpoint[11]=&bv11
    vpoint[12]=&bv12
    vpoint[13]=&bv13
    vpoint[14]=&bv14
    vpoint[15]=&bv15
    vpoint[16]=&bv16
    vpoint[17]=&bv17
    vpoint[18]=&bv18
    vpoint[19]=&bv19
    vpoint[20]=&bv20
    vpoint[21]=&bv21
    vpoint[22]=&bv22
    vpoint[23]=&bv23
    vpoint[24]=&bv24
    vpoint[25]=&bv25
    vpoint[26]=&bv26
    vpoint[27]=&bv27
    vpoint[28]=&bv28
    vpoint[29]=&bv29
    vpoint[30]=&bv30
    vpoint[31]=&bv31
    vpoint[32]=&bv32
    vpoint[33]=&bv33
    vpoint[34]=&bv34
    vpoint[35]=&bv35
    vpoint[36]=&bv36
    vpoint[37]=&bv37
    vpoint[38]=&bv38
    vpoint[39]=&bv39
    vpoint[40]=&bv40
*/  

    laidpreg1=st_matrix("laidpreg1")
    laidpreg2=st_matrix("laidpreg2")
    laidpreg3=st_matrix("laidpreg3")
    laidpreg4=st_matrix("laidpreg4")
    laidpreg5=st_matrix("laidpreg5")
	laidplight1=st_matrix("laidplight1")
    laidplight2=st_matrix("laidplight2")
    laidplight3=st_matrix("laidplight3")
    laidplight4=st_matrix("laidplight4")

    lpreg_1_1=st_matrix("lpreg_1_1")
    lpreg_2_1=st_matrix("lpreg_2_1")
    lpreg_3_1=st_matrix("lpreg_3_1")
    lpreg_4_1=st_matrix("lpreg_4_1")
    lpreg_5_1=st_matrix("lpreg_5_1")
	
	lpreg_1_2=st_matrix("lpreg_1_2")
    lpreg_2_2=st_matrix("lpreg_2_2")
    lpreg_3_2=st_matrix("lpreg_3_2")
    lpreg_4_2=st_matrix("lpreg_4_2")
    lpreg_5_2=st_matrix("lpreg_5_2")
	
	lpreg_1_3=st_matrix("lpreg_1_3")
    lpreg_2_3=st_matrix("lpreg_2_3")
    lpreg_3_3=st_matrix("lpreg_3_3")
    lpreg_4_3=st_matrix("lpreg_4_3")
    lpreg_5_3=st_matrix("lpreg_5_3")
	
	lpreg_1_4=st_matrix("lpreg_1_4")
    lpreg_2_4=st_matrix("lpreg_2_4")
    lpreg_3_4=st_matrix("lpreg_3_4")
    lpreg_4_4=st_matrix("lpreg_4_4")
    lpreg_5_4=st_matrix("lpreg_5_4")
	
	lpreg_1_5=st_matrix("lpreg_1_5")
    lpreg_2_5=st_matrix("lpreg_2_5")
    lpreg_3_5=st_matrix("lpreg_3_5")
    lpreg_4_5=st_matrix("lpreg_4_5")
    lpreg_5_5=st_matrix("lpreg_5_5")
	
	lplight_1_1=st_matrix("lplight_1_1")
    lplight_2_1=st_matrix("lplight_2_1")
    lplight_3_1=st_matrix("lplight_3_1")
    lplight_4_1=st_matrix("lplight_4_1")
	
	lplight_1_2=st_matrix("lplight_1_2")
    lplight_2_2=st_matrix("lplight_2_2")
    lplight_3_2=st_matrix("lplight_3_2")
    lplight_4_2=st_matrix("lplight_4_2")

	lplight_1_3=st_matrix("lplight_1_3")
    lplight_2_3=st_matrix("lplight_2_3")
    lplight_3_3=st_matrix("lplight_3_3")
    lplight_4_3=st_matrix("lplight_4_3")

	lplight_1_4=st_matrix("lplight_1_4")
    lplight_2_4=st_matrix("lplight_2_4")
    lplight_3_4=st_matrix("lplight_3_4")
    lplight_4_4=st_matrix("lplight_4_4")
	
	lp1_1=st_matrix("lp1_1")
	lp1_2=st_matrix("lp1_2")
	lp2_1=st_matrix("lp2_1")
	lp2_2=st_matrix("lp2_2")
	
	lexp1=st_matrix("lexp1")
	lexp2=st_matrix("lexp2")
	indexall=st_matrix("indexall")
	w=st_matrix("w")	
	w_nest=st_matrix("w_nest")	
end

drop _all

mata
    e1j=J(1000,9,0)
    e2j=J(1000,9,0)
    e3j=J(1000,9,0)
    e4j=J(1000,9,0)
    e5j=J(1000,9,0)
    e6j=J(1000,9,0)
    e7j=J(1000,9,0)
    e8j=J(1000,9,0)
    e9j=J(1000,9,0)
	

for(z=1;z<=1000;z++){
        /*makes price coefficient matrix for each draw*/
g=(lpreg_1_1[z,1],lpreg_1_2[z,1],lpreg_1_3[z,1],lpreg_1_4[z,1],lpreg_1_5[z,1]
  \lpreg_2_1[z,1],lpreg_2_2[z,1],lpreg_2_3[z,1],lpreg_2_4[z,1],lpreg_2_5[z,1]  
  \lpreg_3_1[z,1],lpreg_3_3[z,1],lpreg_3_3[z,1],lpreg_3_4[z,1],lpreg_3_5[z,1]  
  \lpreg_4_1[z,1],lpreg_4_4[z,1],lpreg_4_3[z,1],lpreg_4_4[z,1],lpreg_4_5[z,1]  
  \lpreg_5_1[z,1],lpreg_5_5[z,1],lpreg_5_3[z,1],lpreg_5_4[z,1],lpreg_5_5[z,1])

g2=(lplight_1_1[z,1],lplight_1_2[z,1],lplight_1_3[z,1],lplight_1_4[z,1]
   \lplight_2_1[z,1],lplight_2_2[z,1],lplight_2_3[z,1],lplight_2_4[z,1]  
   \lplight_3_1[z,1],lplight_3_3[z,1],lplight_3_3[z,1],lplight_3_4[z,1]  
   \lplight_4_1[z,1],lplight_4_4[z,1],lplight_4_3[z,1],lplight_4_4[z,1])
  
  
/*makes expenditure coefficients*/
b=(laidpreg1[z]\laidpreg2[z]\laidpreg3[z]\laidpreg4[z]\laidpreg5[z])
b2=(laidplight1[z]\laidplight2[z]\laidplight3[z]\laidplight4[z])
delta_mid_reg_reg=(lp1_1[z])
delta_mid_reg_light=(lp1_2[z])
delta_mid_light_reg=(lp2_1[z])
delta_mid_light_light=(lp2_2[z])
delta_top=(indexall[z])
beta_reg=(lexp1[z])
beta_light=(lexp2[z])



sx=J(9,47,0)
	for(j=1;j<=47;j++){           
		sx[.,j]=(*spoint[j])[.,z]
	}	            

	s=mean(sx',1)
	s1=s[1..5]'
	s2=s[6..9]'
	w_all=mean(w) 
	wnest=mean(w_nest)
	ereg=J(5,5,0)
	elight=J(4,4,0)
	ereg_light=J(5,4,0)
	elight_reg=J(4,5,0)

for(m=1;m<=5;m++){
	for(n=1;n<=5;n++){
		ereg[m,n]=(1+delta_mid_reg_reg)*wnest[n]+beta_reg*(1+delta_top)*w_all[n]+g[m,n]/s1[m]+(b[m]/s1[m])*(delta_mid_reg_reg*wnest[n]+w[n]*(1+delta_top)*beta_reg)
	}
}

for(m=1;m<=4;m++){
	for(n=1;n<=4;n++){
		elight[m,n]=(1+delta_mid_light_light)*wnest[n+5]+beta_light*(1+delta_top)*w_all[n+5]+g2[m,n]/s2[m]+(b2[m]/s2[m])*(delta_mid_light_light*wnest[n+5]+w[n+5]*(1+delta_top)*beta_light)
	}
}
ereg=ereg-I(5)
elight=elight-I(4)	



for(m=1;m<=5;m++){
	for(n=1;n<=4;n++){
		ereg_light[m,n]=(1+b[m]/s1[m])*(beta_reg*(1+delta_top)*w[n+5]+delta_mid_reg_light*wnest[n+5])
	}
}

for(m=1;m<=4;m++){
	for(n=1;n<=5;n++){
		elight_reg[m,n]=(1+b2[m]/s2[m])*(beta_light*(1+delta_top)*w[n]+delta_mid_light_reg*wnest[n])
	}
}


		e1j[z,1]=ereg[1,1]
		e1j[z,2]=ereg[1,2]
		e1j[z,3]=ereg[1,3]
		e1j[z,4]=ereg[1,4]
		e1j[z,5]=ereg[1,5]
		e1j[z,6]=ereg_light[1,1]
		e1j[z,7]=ereg_light[1,2]
		e1j[z,8]=ereg_light[1,3]
		e1j[z,9]=ereg_light[1,4]
		
		e2j[z,1]=ereg[2,1]
		e2j[z,2]=ereg[2,2]
		e2j[z,3]=ereg[2,3]
		e2j[z,4]=ereg[2,4]
		e2j[z,5]=ereg[2,5]
		e2j[z,6]=ereg_light[2,1]
		e2j[z,7]=ereg_light[2,2]
		e2j[z,8]=ereg_light[2,3]
		e2j[z,9]=ereg_light[2,4]
		
		e3j[z,1]=ereg[3,1]
		e3j[z,2]=ereg[3,2]
		e3j[z,3]=ereg[3,3]
		e3j[z,4]=ereg[3,4]
		e3j[z,5]=ereg[3,5]
		e3j[z,6]=ereg_light[3,1]
		e3j[z,7]=ereg_light[3,2]
		e3j[z,8]=ereg_light[3,3]
		e3j[z,9]=ereg_light[3,4]
		
		e4j[z,1]=ereg[4,1]
		e4j[z,2]=ereg[4,2]
		e4j[z,3]=ereg[4,3]
		e4j[z,4]=ereg[4,4]
		e4j[z,5]=ereg[4,5]
		e4j[z,6]=ereg_light[4,1]
		e4j[z,7]=ereg_light[4,2]
		e4j[z,8]=ereg_light[4,3]
		e4j[z,9]=ereg_light[4,4]
		
		e5j[z,1]=ereg[5,1]
		e5j[z,2]=ereg[5,2]
		e5j[z,3]=ereg[5,3]
		e5j[z,4]=ereg[5,4]
		e5j[z,5]=ereg[5,5]
		e5j[z,6]=ereg_light[5,1]
		e5j[z,7]=ereg_light[5,2]
		e5j[z,8]=ereg_light[5,3]
		e5j[z,9]=ereg_light[5,4]
		
		e6j[z,1]=elight_reg[1,1]
		e6j[z,2]=elight_reg[1,2]
		e6j[z,3]=elight_reg[1,3]
		e6j[z,4]=elight_reg[1,4]
		e6j[z,5]=elight_reg[1,5]
		e6j[z,6]=elight[1,1]
		e6j[z,7]=elight[1,2]
		e6j[z,8]=elight[1,3]
		e6j[z,9]=elight[1,4]
		
		e7j[z,1]=elight_reg[2,1]
		e7j[z,2]=elight_reg[2,2]
		e7j[z,3]=elight_reg[2,3]
		e7j[z,4]=elight_reg[2,4]
		e7j[z,5]=elight_reg[2,5]
		e7j[z,6]=elight[2,1]
		e7j[z,7]=elight[2,2]
		e7j[z,8]=elight[2,3]
		e7j[z,9]=elight[2,4]
		
		e8j[z,1]=elight_reg[3,1]
		e8j[z,2]=elight_reg[3,2]
		e8j[z,3]=elight_reg[3,3]
		e8j[z,4]=elight_reg[3,4]
		e8j[z,5]=elight_reg[3,5]
		e8j[z,6]=elight[3,1]
		e8j[z,7]=elight[3,2]
		e8j[z,8]=elight[3,3]
		e8j[z,9]=elight[3,4]
		
		e9j[z,1]=elight_reg[4,1]
		e9j[z,2]=elight_reg[4,2]
		e9j[z,3]=elight_reg[4,3]
		e9j[z,4]=elight_reg[4,4]
		e9j[z,5]=elight_reg[4,5]
		e9j[z,6]=elight[4,1]
		e9j[z,7]=elight[4,2]
		e9j[z,8]=elight[4,3]
		e9j[z,9]=elight[4,4]
		
		
}
mean(e1j)
mean(e2j)
mean(e3j)				      
mean(e4j)
mean(e5j)
mean(e6j)
mean(e7j)				      
mean(e8j)
mean(e9j)

aidsolselas=J(9,9,0)

aidsolselas[1,.]=sqrt(mean(e1j:*e1j)-mean(e1j):*mean(e1j))				      
aidsolselas[2,.]=sqrt(mean(e2j:*e2j)-mean(e2j):*mean(e2j))				      
aidsolselas[3,.]=sqrt(mean(e3j:*e3j)-mean(e3j):*mean(e3j))
aidsolselas[4,.]=sqrt(mean(e4j:*e4j)-mean(e4j):*mean(e4j))
aidsolselas[5,.]=sqrt(mean(e5j:*e5j)-mean(e5j):*mean(e5j))				      
aidsolselas[6,.]=sqrt(mean(e6j:*e6j)-mean(e6j):*mean(e6j))				      
aidsolselas[7,.]=sqrt(mean(e7j:*e7j)-mean(e7j):*mean(e7j))				      
aidsolselas[8,.]=sqrt(mean(e8j:*e8j)-mean(e8j):*mean(e8j))				      
aidsolselas[9,.]=sqrt(mean(e9j:*e9j)-mean(e9j):*mean(e9j))				      				      

end





cap log close






mata
    bsfaid=J(1000,9,0)
	regsolution=J(1000,1,0)
    for(z=1;z<=1000;z++){
		g=(lpreg_1_1[z,1],lpreg_1_2[z,1],lpreg_1_3[z,1],lpreg_1_4[z,1],lpreg_1_5[z,1]
		\lpreg_2_1[z,1],lpreg_2_2[z,1],lpreg_2_3[z,1],lpreg_2_4[z,1],lpreg_2_5[z,1]  
		\lpreg_3_1[z,1],lpreg_3_3[z,1],lpreg_3_3[z,1],lpreg_3_4[z,1],lpreg_3_5[z,1]  
		\lpreg_4_1[z,1],lpreg_4_4[z,1],lpreg_4_3[z,1],lpreg_4_4[z,1],lpreg_4_5[z,1]  
		\lpreg_5_1[z,1],lpreg_5_5[z,1],lpreg_5_3[z,1],lpreg_5_4[z,1],lpreg_5_5[z,1])

		g2=(lplight_1_1[z,1],lplight_1_2[z,1],lplight_1_3[z,1],lplight_1_4[z,1]
		\lplight_2_1[z,1],lplight_2_2[z,1],lplight_2_3[z,1],lplight_2_4[z,1]  
		\lplight_3_1[z,1],lplight_3_3[z,1],lplight_3_3[z,1],lplight_3_4[z,1]  
		\lplight_4_1[z,1],lplight_4_4[z,1],lplight_4_3[z,1],lplight_4_4[z,1])
  
  
		/*makes expenditure coefficients*/
		b=(laidpreg1[z]\laidpreg2[z]\laidpreg3[z]\laidpreg4[z]\laidpreg5[z])
		b2=(laidplight1[z]\laidplight2[z]\laidplight3[z]\laidplight4[z])
		delta_mid_reg_reg=(lp1_1[z])
		delta_mid_reg_light=(lp1_2[z])
		delta_mid_light_reg=(lp2_1[z])
		delta_mid_light_light=(lp2_2[z])
		delta_top=(indexall[z])
		beta_reg=(lexp1[z])
		beta_light=(lexp2[z])
		delta=(delta_mid_reg_reg, delta_mid_reg_light\delta_mid_light_reg, delta_mid_light_light)
		beta=(beta_reg\beta_light)

		temp=J(47,9,0)
		check=J(47,1,0)	
        for(j=1;j<=47;j++){
            p=(*ppoint[j])[.,z]
            pinit=(*ppoint[j])[.,z]
            s=(*spoint[j])[.,z]
            aidomega=J(9,9,0)
		
			for(x=1;x<=5;x++){
				aidomega[x,x]=(-1+(1+delta_mid_reg_reg)*w_nest[j,x]+beta_reg*(1+delta_top)*w[j,x]+g[x,x]/s[x]+(b[x]/s[x])*(delta_mid_reg_reg*w_nest[j,x]+w[j,x]*(1+delta_top)*beta_reg))*X[j,1]*s[x]
			}
			for(x=6;x<=9;x++){
				aidomega[x,x]=(-1+(1+delta_mid_light_light)*w_nest[j,x]+beta_light*(1+delta_top)*w[j,x]+g2[x-5,x-5]/s[x]+(b2[x-5]/s[x])*(delta_mid_light_light*w_nest[j,x]+w[j,x]*(1+delta_top)*beta_light))*X[j,2]*s[x]
			}
			aidomega[6,1]=((1+b[1]/s[1])*(beta_reg*(1+delta_top)*w[j,6]+delta_mid_reg_light*w_nest[j,6]))*X[j,1]*s[1]			            			
			aidomega[1,6]=((1+b2[1]/s[6])*(beta_light*(1+delta_top)*w[j,1]+delta_mid_light_reg*w_nest[j,1]))*X[j,2]*s[6]
			aidomega[7,2]=((1+b[2]/s[2])*(beta_reg*(1+delta_top)*w[j,7]+delta_mid_reg_light*w_nest[j,7]))*X[j,1]*s[2]
			aidomega[2,7]=((1+b2[2]/s[7])*(beta_light*(1+delta_top)*w[j,2]+delta_mid_light_reg*w_nest[j,2]))*X[j,2]*s[7]
			aidomega[8,3]=((1+b[3]/s[3])*(beta_reg*(1+delta_top)*w[j,8]+delta_mid_reg_light*w_nest[j,8]))*X[j,1]*s[3]
			aidomega[3,8]=((1+b2[3]/s[8])*(beta_light*(1+delta_top)*w[j,3]+delta_mid_light_reg*w_nest[j,3]))*X[j,2]*s[8]
			aidomega[9,4]=((1+b[4]/s[4])*(beta_reg*(1+delta_top)*w[j,9]+delta_mid_reg_light*w_nest[j,9]))*X[j,1]*s[4]
			aidomega[4,9]=((1+b2[4]/s[9])*(beta_light*(1+delta_top)*w[j,4]+delta_mid_light_reg*w_nest[j,4]))*X[j,2]*s[9]
			
			Xi=J(9,9,0)
			Xi[1,1]=X[j,1]
			Xi[2,2]=X[j,1]
			Xi[3,3]=X[j,1]
			Xi[4,4]=X[j,1]
			Xi[5,5]=X[j,1]
			Xi[6,6]=X[j,2]
			Xi[7,7]=X[j,2]
			Xi[8,8]=X[j,2]
			Xi[9,9]=X[j,2]

			mu=-luinv(aidomega)*Xi*s
			aidcost=p:*(J(9,1,1)-mu)
			pinit=(*ppoint[j])[.,z]
		
			a=csolve(&msyraids(),pinit,1e-6,400,s,p,aidcost,g,g2,b,b2,w[j,.]',w_nest[j,.]',X[j,.]',delta,beta,delta_top)		
			ws=msyraids(a[1..9,.],s,p,aidcost,g,g2,b,b2,w[j,.]',w_nest[j,.]',X[j,.]',delta,beta,delta_top)
			check[j,1]=((a[10,1]+(a[1..rows(a)-1,.]==J(rows(a)-1,1,.))+(ws[1,.]==.)+(ws[2,.]==.)+(ws[3,.]==.)+(ws[4,.]==.)+(ws[5,.]==.)+(ws[6,.]==.)+(ws[7,.]==.)+(ws[8,.]==.)+(ws[9,.]==.))==0)			
			a=a[1..9,.]
			a=(a-p):/p
			temp[j,.]=a'
    }
			bsfaid[z,.]=mm_median(select(temp,check))    		      
}       
	st_matrix("bsfaid",bsfaid)
end

svmat bsfaid
log using syrtable4col4, replace
	centile bsfaid1, centile(2.5, 97.5)
	centile bsfaid2, centile(2.5, 97.5)
	centile bsfaid3, centile(2.5, 97.5)
	centile bsfaid4, centile(2.5, 97.5)
	centile bsfaid5, centile(2.5, 97.5)
	centile bsfaid6, centile(2.5, 97.5)
	centile bsfaid7, centile(2.5, 97.5)
	centile bsfaid8, centile(2.5, 97.5)
	centile bsfaid9, centile(2.5, 97.5)
	log close
	
save c:/data/restat_12_9/bootstrap_results/post_bsivaids, replace
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


forval x=1/9{
	replace bsfaid`x'=100*bsfaid`x'
}

rename bsfaid1 Aunt_Jemima
rename bsfaid2 Hungry_Jack
rename bsfaid3 Log_Cabin
rename bsfaid4 Mrs_Butterworth
rename bsfaid5 PrivateLabel

xxxxxxxxxxxxxxxxxxx
twoway histogram Castrol, title("Castrol") xtitle("Percent Price Change") saving(hc, replace)
twoway hist Havoline, title("Havoline") xtitle("Percent Price Change") saving(hh, replace)
twoway histogram Mobil, title("Mobil") xtitle("Percent Price Change") saving(hm, replace)
twoway histogram Pennzoil, title("Pennzoil") xtitle("Percent Price Change") saving(hpe, replace)
twoway hist PrivateLabel, title("Private Label") xtitle("Percent Price Change") saving(hp, replace)
twoway hist Quaker, title("Quaker") xtitle("Percent Price Change") saving(hq, replace)
twoway hist Valvoline, title("Valvoline") xtitle("Percent Price Change") saving(hv, replace)

graph combine hc.gph hh.gph hm.gph hpe.gph hp.gph hq.gph hv.gph, title(Sampling Distributions of Price Changes)
graph export olsaidsoil.eps, replace

/*Are sampling distributions normal?*/
/*qq plots*/
qnorm Castrol, title("Castrol") saving(qc, replace)
qnorm Havoline, title("Havoline") saving(qh, replace)
qnorm Mobil, title("Mobil") saving(qm, replace)
qnorm Pennzoil, title("Pennzoil") saving(qp, replace)
qnorm PrivateLabel, title("Private Label") saving(qplo, replace)
qnorm Quaker, title("Quaker State") saving(qq, replace)
qnorm Valvoline, title("Valvoline") saving(qv, replace)


graph combine qc.gph qh.gph qm.gph qp.gph qplo.gph qq.gph qv.gph
graph export qolsaidsoil.eps, replace

/*Jarque-Bera Tests for Normality*/
egen sc=skew(Castrol)
egen kc=kurt(Castrol)
gen jbc=_N*(saj^2/6+(kaj-3)^2/24)

egen sh=skew(Havoline)
egen kh=kurt(Havoline)
gen jbh=_N*(sh^2/6+(kh-3)^2/24)

egen sm=skew(Mobil)
egen km=kurt(Mobil)
gen jbm=_N*(sm^2/6+(km-3)^2/24)

egen sp=skew(Pennzoil)
egen kp=kurt(Pennzoil)
gen jbp=_N*(sp^2/6+(kp-3)^2/24)

egen splo=skew(PrivateLabel)
egen kplo=kurt(PrivateLabel)
gen jbplo=_N*(splo^2/6+(kplo-3)^2/24)

egen sq=skew(QuakerState)
egen kq=kurt(QuakerState)
gen jbq=_N*(sq^2/6+(kq-3)^2/24)

egen sv=skew(Valvoline)
egen kv=kurt(Valvoline)
gen jbv=_N*(sv^2/6+(kv-3)^2/24)

sum jbaj jbhj jblc jbmb jbplo
save oilaidsols, replace









xxxxx
